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to  correlate  the  optimal  sizes  of  bodies  with  specified  external  forced  convection  heat  transfer, 
when  the  objective  is  to  minimize  the  total  rate  of  entropy  generation,  or  the  total  cost. 
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Abstract 

This  research  project  focused  on  the  friction  and  heat  transfer  characteristics  of  a  solid 
surface  covered  with  flexible  fibers.  The  work  consisted  of  analysis,  numerical  simulaitons,  and 
laboratory  experiments.  Related  developments  that  grew  out  of  this  work  were  investigated  in 
the  later  part  of  the  project:  specifically,  the  optimal  spacings  between  cylinders  (e.g.  fibers)  with 
heat  transfer,  and  the  optimal  sizes  of  bodies  immersed  in  an  external  flow. 

Chapter  1  documents  the  pressure  drop  and  heat  transfer  through  bundles  of  parallel 
cylinders  in  a  domain  that  has  been  overlooked,  namely,  low  Reynolds  numbers,  arrays  that  are 
long  in  the  direction  of  flow  such  that  the  flow  is  hydraulically  and  thermally  fully  developed, 
and  cylinders  inclined  relative  to  the  flow  direction.  The  numerical  results  cover  the  range 
l<ReD<30,  0.72<Pr<100,  0.6<(t)<0.95  and  0°  <(3<60%  where  ^  is  the  porosity  of  the 
bundle  as  a  saturated  porous  medium,  and  P  is  the  angle  between  the  cylinder  centerline  and  the 
direction  perpendicular  to  the  flow  direction.  The  accuracy  is  verified  by  means  of  experimental 
measurements  of  the  pressure  drop  across  a  bundle  with  115  cylinders  in  the  flow  direction,  in 
the  range  8<ReD<50, 0.84 <(!>< 0.92  and  0°<P<60°.  The  results  show  that  significant  errors 
may  occur  if  the  available  large-Reo  information  is  extrapolated  to  the  domain  covered  by  this 
study. 

The  heat  transfer,  friction  and  mechanical  (elastic)  interaction  between  an  external 
laminar  flow  and  a  solid  surface  covered  by  a  layer  of  fibers  is  documented  numerically  in 
Chapter  2.  The  flow  is  initially  perpendicular  to  the  surface,  and  the  fibers  can  bend.  In  the  first 
part  of  the  chapter,  it  is  assumed  that  the  fibers  are  inflexible.  It  is  shown  that  the  effect  of  the 
fiber  layer  on  the  overall  heat  transfer  can  be  correlated  in  terms  of  the  fraction  of  the  external 
flow  that  penetrates  into  the  fiber  layer.  The  second  part  focuses  on  the  effect  of  fiber  bending, 
which  is  described  by  a  new  dimensionless  group:  the  stiffness  number  S  =  El  /  (pU^L^).  It  is 

shown  that  the  wall  heat  transfer,  friction  and  fiber  layer  flow  fraction  exhibit  a  sudden  decrease 
when  S  drops  below  a  critical  value.  Sc-  The  critical  stiffness  number  can  be  correlated  as 
Sc  =  C(H/L)5  (D/H)2  /  (1  -  $)  where  C  is  a  constant  of  order  0.4,  H  and  D  are  the  fiber  length 
and  diameter,  L  is  the  half-length  of  the  solid  wall,  and  (])  is  the  porosity  of  the  fiber  layer. 

Chapter  3  summarizes  a  theoretical,  numerical  and  experimental  study  of  how  to  select 
the  spacing  (S)  between  horizontal  cylinders  in  an  array  with  laminar  natural  convection,  such 
that  the  total  heat  transfer  (q)  between  the  array  and  the  ambient  is  maximized.  The  volume 
occupied  by  the  array  (height  H,  width  W,  cylinder  length  L)  and  the  cylinder  diameter  (D)  are 
arbitrary  but  fixed,  while  the  spacing  (or  number  of  cylinders  in  the  array)  varies.  The  optimal 
spacing  and  maximum  heat  transfer  results  predicted  theoretically  are  developed  into  accurate 
and  well  tested  correlations  by  means  of  numerical  simulations  and  experimental  measurements. 
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The  recommended  correlations  are  S^pj/D  =  2.72  +  0.263  and 

q^ax  =  0.448[(H/D)^''^RaD^'^  where  q^ax  is  the  dimensionless  maximum  overall  thermal 

conductance,  q^^ax  =  Tmax  LW  k(Tw  -  Too)].  The  optimal  spacing  is  relatively  insensitive 
to  whether  the  cylinders  are  isothermal  or  with  uniform  heat  flux.  The  corresponding  problem  of  k 

finding  the  optimal  spacing  when  the  heat  transfer  is  by  forced  convection  is  considered  in 
Chapter  4. 

Chapter  5  shows  how  to  correlate  the  optimal  sizes  of  bodies  with  specified  external 
forced  convection  heat  transfer,  when  the  objective  is  to  minimize  (1)  the  total  rate  of  entropy 
generation,  or  (2)  the  total  cost.  Specific  results  are  reported  for  the  cylinder  in  cross-flow,  the 
sphere,  and  the  flat  plate  in  parallel  flow.  The  results  are  presented  in  terms  of  dimensionless 
parameters  that  show  that  regardless  of  the  body  shape,  the  optimal  size  increases  monotonically 
with  the  heat  transfer  duty  parameter.  A  universal  correlation  is  obtained  by 
nondimensionalizing  the  body  size  as  a  Reynolds  number  based  on  the  transversal  length  scale  of 
the  flow,  and  by  replacing  the  body  with  an  equivalent  one  in  two-dimensional  flow.  The 
optimal  sizes  correlated  in  this  manner  can  be  applied  to  bodies  of  other  shapes. 
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Chapter  1 


FORCED  CONVECTION  IN  BANKS  OF  INCLINED  CYLINDERS 
AT  LOW  REYNOLDS  NUMBERS 


1.1  Introduction 

This  chapter  documents  the  flow  and  heat  transfer  characteristics  of  bundles  of  parallel 
cylinders  that  are  inclined  relative  to  the  free  stream  that  bathes  them.  The  emphasis  is  on  (a)  the 
low  Reynolds  number  range  Reo  ^  30,  (b)  the  fully  developed  regime,  and  (c)  the  effect  of  the 
angle  of  inclination  on  the  pressure  drop  and  heat  transfer  coefficient.  These  aspects  have  been 
overlooked  in  spite  of  the  large  research  effort  that  has  been  devoted  to  banks  of  cylinders  in 
cross-flow.  The  progress  on  this  general  topic  has  been  reviewed  by  Zukauskas  (1972, 1987a,  b) 
and  Kays  and  London  (1984). 

In  his  more  recent  review,  Zukauskas  (1987b)  proposed  a  correction  factor  to  account  for 
the  effect  of  angle  of  inclination  on  pressure  drop  calculations.  It  is  unclear  whether  his  method 
is  universally  valid,  or  applies  only  to  a  certain  range  of  Reynolds  numbers  and  a  certain  type  of 
cylinder  array.  One  objective  of  the  present  study  was  to  clarify  the  conditions  under  which 
Zukauskas’s  correction  is  applicable. 

The  flow  and  heat  transfer  results  assembled  in  this  chapter  were  developed  in  three 
phases.  First,  numerous  cases  of  the  flow  and  temperature  fields  in  long  bundles  of  inclined 
cylinders  were  simulated  numerically,  based  on  detailed  three-dimensional  calculations.  In  the 
second  phase,  the  heat  and  fluid  flow  results  were  correlated  into  a  saturated  porous  medium 
model  that  accounts  for  the  angle  of  inclination  and  the  low  Reynolds  number  range.  The  third 
phase  was  experimental,  where  the  accuracy  of  the  numerical  results  was  tested  against 
laboratory  measurements  of  the  pressure  drops  through  bundles  at  several  cylinder  inclinations 
and  Reynolds  numbers. 

1.2  Mathematical  Formulation 

The  numerical  part  of  the  study  refers  to  an  array  of  staggered  cylinders,  with  their 
centers  arranged  in  an  array  of  equilateral  triangles.  Two  orientations  of  the  flow  relative  to  the 
array  were  examined  (Fig.  1.1,  top).  The  entire  array  was  tilted  to  an  angle  p  that  varied  from  0° 
to  60°  relative  to  the  perpendicular  (Fig.  1.1,  bottom).  In  the  numerical  simulations  (and  unlike 
in  the  experiments,  section  1.6)  the  spacing  between  cylinders  was  fixed  so  that  the  porosity  of 
the  bundle  ^  was  not  a  function  of  the  angle  p.  In  an  equilateral  triangle  array  the  porosity  is  <{>  = 
1  -  0.907  (D/S)2,  where  S  is  the  distance  between  cylinder  centers. 
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Figure  1.1  Arrays  of  parallel  cylinders  with  forced  convection  heat  transfer  (top),  and  three- 
dimensional  computational  domain  (bottom). 
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In  the  range  Reo  ^  30  the  flow  is  laminar  and  the  wake  behind  each  cylinder  is  steady  and 
symmetric.  This  allowed  us  to  view  the  array  as  a  sandwich  of  many  channels  like  the  one 
highlighted  in  the  middle  of  Fig.  1.1,  and  illustrated  by  the  dashed  lines  in  the  top  drawings  of 
Fig.  1.1.  There  was  no  mixing  between  the  channels.  The  computational  domain  contained  only  one 
channel. 

The  equations  that  account  for  the  conservation  of  mass,  momentum  and  energy  in  the 
fluid  regions  of  the  three-dimensional  frame  defined  in  Fig.  1.1  are 


3u  I  dv 
dx^dy 


(1.1) 


u.^  =  _i3P+v^ 

“)  3xj  P  ''  3x  J 


(1.2) 


dx  dy  dz 


(1.3) 


where  (u,v,w)  are  the  velocity  components  aligned  with  the  (x,y,z)  axes.  Equation  (1.2)  makes 
use  of  the  tensor  notation  convention  that  triple  summation  is  indicated  by  repeated  indices.  The 
boundary  conditions  on  the  cylinder  surfaces  account  for  no  slip,  no  penetration  and  uniform 
temperature. 


Ui  =  0  and  T  =  Tq  on  all  cylinder  surfaces  (1.4) 

Along  the  top  and  bottom  of  the  channel,  and  on  the  side-wall  regions  contained  between 
cylinders  the  conditions  were  zero  shear,  no  penetration  and  zero  heat  flux. 


3u 

dz 

=  0, 

^  =  0 
dz 

W  = 

=  0 

and 

II 

o 

at  z  =  0,H 

(1.5) 

QJ\CV 

=0, 

v  =  0. 

II 

0, 

and 

o 

II 

at  y  =  0,W 

(1.6) 

The  flow  was  isothermal  and  longitudinally  uniform  at  the  channel  inlet. 


u  =  Uoo,  v  =  0,  w  =  0  and  T  =  Too  at  x  =  0  (1.7) 

The  channel  outlet  conditions  were  zero  stress  and  zero  longitudinal  heat  flux, 

’‘"L  (1.8) 

The  equations  and  boundary  conditions  were  nondimensionalized  by  using  D,  Uoo  and  (Tq  -  Too) 
as  representative  scales. 


(X,Y,Z)  = 


(x,y,z) 
D  ’ 


(U,V,W)  = 


(u,v,w) 


U 


oo 


(1.9) 
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(1.10) 


T-T 

0  =  -  “ 


P 


To-Too  pU^ 

The  nondimensional  counterparts  of  equations  (1.1)-(1.8),  which  are  omitted  for  brevity,  contain 
the  nondimensional  numbers 

The  system  of  nondimensional  equations  was  solved  using  the  finite  element  package 
FIDAP  (1991).  When  the  cylinders  are  inclined  (P  ^  0)  the  problem  is  inherently  three¬ 


dimensional  and  thus  a  3-D  solver  was  required.  The  computational  domain  was  built  up  of  8 
node  brick  elements.  A  Stokes  flow  solution  was  used  to  initialize  the  solution  for  the  mass  and 
momentum  equations  when  P  =  0°,  and  then  the  solution  to  each  prior  run  was  used  to  initialize 
the  next  ran,  as  P  was  increased.  Since  the  energy  equation  is  only  weakly  coupled  to  the  mass 
and  momentum  equations,  it  was  solved  separately  for  the  various  Pr  values  after  the  velocity 
field  was  determined. 

Numerical  testing  indicated  that  a  channel  with  40  cylinders  in  the  longitudinal  direction 
allowed  for  pressure  measurements  at  rows  11,  19  and  27  that  were  not  effected  by  the  entrance 
or  exit  regions.  We  verified  that  the  flow  was  fully  developed  by  calculating  the  pressures  at 
three  equidistant  longitudinal  locations  (say  a,  b  and  c),  and  checking  that  the  pressure  drops 
matched.  Pa  -  Pb  =  Pb  -  Pc-  In  all  cases,  the  discrepancy  between  the  two  pressure  drops  was 
less  than  1  percent. 

The  effect  of  doubling  the  channel  height  H  was  found  to  be  insignificant  (less  than 
1  percent  change  in  the  pressure  gradient)  if  the  H  >  6  and  a  free  slip  boundary  condition  was 
used  at  the  top  and  bottom  of  the  channel.  The  solution  proved  to  be  extremely  sensitive  to 
changes  in  channel  height  when  a  no  slip  condition  was  imposed  on  the  top  and  bottom  planes  of 
the  channel.  Convergence  testing  for  the  mesh  density  was  performed  with  p  =  45°.  The  number 
of  elements  was  doubled  in  each  direction  until  such  doubling  produced  a  change  in  the  pressure 
gradient  of  less  than  3  percent.  These  tests  were  performed  for  each  porosity,  and  for  both 
orientations  of  the  array. 

All  of  the  final  numerical  meshes  for  determining  the  pressure  gradient  required  between 
20,000  and  35,000  elements.  Typically  the  initial  ran,  which  requires  a  Stokes  flow  solution  for 
initialization,  required  10  to  20  minutes  of  CPU  time  on  a  Cray  Y-MP.  The  runs  that  followed 
generally  required  about  10  minutes  of  CPU  time  each.  The  temporary  space  required  for  the 
matrix  solver  ranged  from  500  to  several  thousands  of  megabytes.  This,  in  combination  with  the 
increased  computational  cost  of  doing  proper  convergence  testing,  made  a  supercomputer 
necessary.  The  convergence  testing  for  stability  in  the  calculated  Nusselt  number  indicated  that 
twice  as  many  elements  were  necessary  in  the  x  direction  to  achieve  a  3  percent  stability  in  the 
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Nusselt  number  than  were  needed  for  the  same  stability  in  the  pressure  gradient.  To  solve  for  the 
velocity  field  on  the  finer  mesh  required  about  50  minutes  of  CPU  time  on  the  Cray  Y-MP. 
Three  additional  minutes  were  needed  to  solve  the  heat  transfer  problem  for  each  Prandtl 
number. 

The  Nusselt  number  calculation  could  not  be  performed  for  the  entire  range  defined  by 
0.6  <  (()  <  0.95,  0.72  <  Pr  <  100  and  1  <  Reo  ^  30.  For  example,  in  the  case  4)  =  0.6,  Pr  =  0.72, 
and  ReD  =  10,  the  fluid  reached  thermal  equilibrium  with  the  solid  almost  immediately,  i.e.  after 
the  first  row.  In  the  opposite  extreme  (e.g.  <{)  =  0.95,  Pr  =  100,  ReD  =  10)  the  flow  was  not  fully 
developed  thermally  at  the  downstream  end  of  the  computational  domain,  and  no  fully  developed 
Nu  could  be  reported.  As  a  rule,  if  Nu  varied  by  more  than  5  percent  over  the  20  inner  rows  of 
cylinders,  the  thermal  development  was  considered  incomplete. 

1.3.  Pressure  Drop  Results 

The  results  that  document  the  effect  of  the  cylinder  angle  of  inclination  on  pressure  drop 
are  presented  in  Figs.  1.2a-d.  The  ordinate  shows  the  pressure  gradient  in  an  array  with  inclined 
cylinders,  as  a  fraction  of  the  pressure  gradient  when  the  same  array  is  perpendicular  to  the  same 
flow.  The  first  three  frames  show  that  when  Reo  ^  10  the  relative  pressure  drop  AP(P)/AP(0°)  is 
influenced  only  by  the  angle  P,  and  not  by  the  other  parameters  (Reo,  4),  orientation).  The  P 
effect  can  be  significant;  for  example,  in  all  four  frames  of  Fig.  1.2  the  pressure  drop  decreases  to 
60  percent  of  its  original  value  as  P  increases  from  0'  to  60°. 

Figures  1.2a  and  1.2d  indicate  that  at  Reynolds  numbers  larger  than  10  the  relative 
pressure  drop  depends  on  Rod  and  4),  in  addition  to  p.  The  pressure  drop  is  considerably  more 
sensitive  to  changes  in  the  angle  P  at  low  porosities  (4)  =  0.6)  than  at  high  porosities  (4)  =  0.9). 
When  the  porosity  is  high  the  angle  of  inclination  has  no  visible  effect  on  the  pressure  drop  if  P 
<30°. 

These  results  can  be  compared  with  the  single  relative  pressure  drop  curve  reported  by 
Zukauskas  (1987),  which  falls  right  on  top  of  the  curve  drawn  in  Fig.  2a  for  Reo  =  30  and 
4)  =  0.6.  The  present  results  suggest  that  the  AP(P)/AP(0°)  curve  is  not  universal,  i.e.  it  is  not 
independent  of  Rep  and  4).  This  lack  of  universality  is  particularly  evident  in  the  direction  of 
increasing  Reo,  which  is  the  domain  of  heat  exchanger  applications.  At  the  same  time,  the 
present  results  show  that  a  universal  relative  pressure  drop  curve  different  than  Zukauskas's 
exists  when  Rod  is  smaller  than  10. 

We  also  correlated  the  pressure  drop  data  for  cross-flow  (P  =  0°)  by  modelling  the  array 
and  the  fluid  as  a  saturated  porous  medium.  In  the  small-ReD  limit  the  flow  is  expected  to  follow 
the  Darcy  model. 
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Figure  1.2  The  effect  of  |3,  <|),  Reo  and  array  orientation  on  the  longitudinal  pressure  drop  in  fully  developed  flow. 
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The  permeability  K  was  determined  by  fitting  the  numerical  pressure  gradient  results  to  equation 
(1.12),  and  correlating  the  K  values  using  a  formula  similar  to  the  Carman-Kozeny  relation 
(Carman,  1937), 

125[(i_(,,)2 

with  r  =  0.85  for  orientation  1,  and  r  =  0.8  for  orientation  2.  Equation  (1.13)  agrees  within  6 
percent  with  the  numerical  K  values  deduced  from  equation  (1.12)  for  orientation  1  and  ^  <  0.9. 
The  agreement  between  equations  (1.13)  and  (1.12)  is  within  35  percent  for  all  the  runs  made  for 
orientation  2:  The  average  error  for  all  the  runs  in  both  orientations  was  10  percent.  The 
agreement  for  the  runs  in  orientation  2  is  illustrated  in  Fig.  1.3,  where  the  numbers  listed  after 
each  point  represent  Reo  and  (]).  The  use  of  Kl72  as  length  scale  in  the  dimensionless  groups  on 
the  ordinate  and  abscissa  of  Fig.  1.3  is  standard  in  the  field  of  convection  in  porous  media 
(Cheng,  1978;  Nield  and  Bejan,  1992).  Equation  (1.12)  is  not  shown  in  the  figure:  it  would  be 
represented  by  a  line  of  slope  -1  passing  through  the  point  (1,1),  i.e.  a  line  that  would  pass  right 
through  the  shown  data. 


r 

(1.13) 


1.4  Heat  Transfer  Results  for  Cylinders  in  Cross-Flow 

The  results  of  the  heat  transfer  part  of  the  problem  defined  in  section  1.2  are  summarized 
as  an  average  Nusselt  number 

Nu  =  i^  (1.14) 

The  rate  of  heat  transfer  between  the  fluid  and  the  wetted  area  A  in  any  control  volume  of  the 
computational  domain  is 

q  =  h  A  ATim  =  ih  Cp  (Tout~  Tjj,)  (1.15) 

where  rh  is  the  mass  flowrate  through  the  domain  (Fig.  1.1,  bottom).  The  relation  between  Nu 
and  the  nondimensional  formulation  presented  in  section  1.2  is 


where  A  =  A/D^  and 
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(1.16) 
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Figure  1.3  Porous-medium  correlation  of  the  pressure  drop  results  for  crossflow. 
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Experimental  studies  usually  report  Nu  averaged  over  the  entire  array,  because  that  value  is 
experimentally  accessible.  We  calculated  the  "local"  Nusselt  number  associated  with  (averaged 
over)  a  certain  row  of  cylinders  in  order  to  monitor  the  thermal  development  of  the  flow,  from 
one  row  to  the  next. 

The  results  obtained  for  P  =  0°  show  that  in  general  the  row  Nusselt  number  Nu  in  the 
first  few  rows  is  larger  than  in  later  rows.  This  trend  is  illustrated  in  Fig.  1.4.  We  found  that  the 
Nu  value  for  the  first  row  was  very  sensitive  to  grid  refinement.  After  a  few  rows  Nu 
approached  a  constant,  which  was  less  sensitive  to  grid  refinement.  Our  convergence  tests  were 
based  on  the  stability  of  the  constant  Nu  value  found  between  rows  7  through  31. 

The  only  correlation  for  low  Reynolds  number  heat  transfer  for  staggered  tube  banks  in 
cross-flow  appears  to  be  a  formula  due  to  Zukauskas  (1987a): 

Nu  =  1.04  Re?)j (1.19) 

where  Nu  is  the  value  for  "inner  rows",  and  ReD,f  is  based  on  the  velocity  averaged  over  the 
minimum  cross  section.  Equation  (1.19)  was  recommended  for  the  range  1.6  <  Reo.f  <  40.  The 
subscripts  f  and  w  indicate  that  the  properties  must  be  evaluated  at  the  flow  mean  temperature 
and,  respectively,  the  cylinder  surface  temperature. 

Our  results  for  fully  developed  cross-flow  in  the  range  1  <  Re^  ^  30  and  0.72  <  Pr  <  100 
are  correlated  by  the  power  law 

Nu  =  c  (Reo  Pr)"  (1.20) 

for  which  the  constants  c  and  n  are  reported  in  Table  1.1.  The  agreement  between  the  numerical 
Nu  values  and  the  values  obtained  based  on  equation  (20)  is  within  4  percent. 

The  most  important  difference  between  the  correlations  (1.19)  and  (1.20),  is  the  smaller 
exponent  n  (i.e.  weaker  Re©  and  Pr  effect)  present  in  equation  (1.20).  To  understand  these 
differences  we  examined  the  experimental  and  numerical  work  on  which  equation  (1.19)  is  based 
(Omohundro  et  al.,  1949;  Bergelin  et  al.,  1949,  1950).  For  example,  equation  (1.20) 
underestimates  by  a  factor  of  2  the  Nusselt  number  found  experimentally  by  Bergelin  et  al. 
(1950)  for  Pr  s  500.  The  numerical  work  of  Chang  et  al.  (1989),  which  was  performed  for 
Pr  =  0.7,  agrees  with  the  Bergelin  et  al.  experiments  when  a  PrO.33  correction  factor  is  applied  to 
the  Nusselt  number.  The  agreement  between  Bergelin  et  al.  and  Chang  et  al.  is  due  to  the 
shortness  of  the  arrays  studied  (10  rows  for  Bergelin  et  al.,  and  5  rows  for  Chang  et  al.),  and 
explains  why  equations  (1.19)  and  (1.20)  are  different. 
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Table  1.1.  Constants  for  the  fully  developed  Nusselt  number  correlation  (1.20). 


<1) 

c 

n 

0.6 

5.64 

0.14 

0.9 

2.34 

0.18 

0.95 

2.0 

0.18 

Figure  1.4  illustrates  the  effect  of  Prandtl  number  on  the  local  Nusselt  number  over  the 
first  30  rows.  When  we  examine  our  numerical  results  for  the  first  ten  cylinder  rows  we  find  that 
the  Nusselt  number  averaged  over  ten  rows  scales  as  PrO-3,  and  that  our  numerical  simulation 
agrees  with  the  findings  of  Bergelin  et  al.  to  within  10  percent.  Nonetheless  at  later  rows  where 
the  flow  is  fully  developed,  our  correlation  (20)  applies.  For  flow  systems  in  which  the  flow 
becomes  fully  developed  early  in  the  array,  it  is  clear  that  Zukauskas's  formula  (1.19)  would 
result  in  significant  error  (although  it  is  correct  for  short  arrays  in  which  entrance  effects 
dominate). 

The  same  explanation  (fully  developed  flow  vs.  entrance  flow)  holds  for  the  different 
exponents  on  the  Reynolds  number  in  equations  (1.19)  and  (1.20).  That  the  Nusselt  number  in 
the  fully  developed  region  is  proportional  to  Reo"  where  n  is  significantly  smaller  than  in 
Zukauskas's  correlation  (1.19)  is  further  supported  by  the  experiments  of  Minakami  et  al.  (1993). 
These  authors  tested  an  array  of  in-line  square  pin  fins  with  10  to  20  rows  in  the  longitudinal 
direction,  and  found  that  Nu  ~  Reo",  where  n  is  between  0.2  and  0.25  in  the  laminar  regime  and 
for  Pr  =  7  (water).  Comparison  of  our  findings  with  the  results  of  Minakami  et  al.  is  reasonable 
because  the  shape  of  the  pins  (square  or  cylindrical)  should  have  only  a  small  effect  in  the  low 
Rod  regime.  Furthermore,  Zukauskas  (1987a)  indicates  that  Nu  values  for  staggered  and  in-line 
arrays  have  the  same  ReD  dependence  for  Reo  >  40. 

1.5.  Heat  Transfer  Results  for  Inclined  Cylinders 

Zukauskas  (1972)  proposed  the  single  curve  of  Fig.  1.5  (top),  as  a  heat  transfer  method  of 
accounting  for  the  effect  of  the  angle  of  inclination  j3  between  the  cylinder  axis  and  the  direction 
perpendicular  to  the  flow.  He  did  not  indicate  the  ReD  range  and  array  dimensions  for  which  this 
curve  is  valid.  Our  results  for  fully  developed  flow,  1  <  Reo  ^  30,  0.72  <  Pr  <  100,  and 
orientation  1  indicate  that  the  angle  effect  on  Nu  is  considerably  smaller  than  indicated  by 
Zukauskas. 

Groehn  (1981)  proposed  the  principle  of  independence  as  a  method  of  accounting  for  the 
effect  of  P  on  the  Nusselt  number.  In  this  method  the  relevant  velocity  scale  is  that  normal  to  the 
cylinder,  hence  the  Reynolds  number  that  should  be  used  in  equation  (1.20)  is  Reo  cos  p.  If  the 
principle  of  independence  is  valid  then  the  ratio  Nu(P)  /  Nu(0°)  should  be  a  function  of  (|)  and  P 


11 


only.  Figure  1.5  (bottom)  illustrates  the  curves  predicted  by  the  principle  of  independence  in 
combination  with  equation  (1.20). 

Figures  1.6a-c  show  how  the  ratio  Nu(p)/Nu(0°)  responds  to  changes  in  P,  (1>  and  Pr  when 
Reo  is  fixedAn  interesting  feature  is  the  evolution  of  the  Pr  =  7  curve  (e.g.  water)  as  the 
porosity  increases  from  0.6  to  0.95.  The  Nu(p)/Nu(0°)  curve  moves  upward  as  ({>  increases,  and 

this  means  that  the  angle  effect  is  weaker  in  a  sparse  array  than  in  a  dense  array.  In  particular, 
when  ([)  =  0.95  and  0°  <  P  <  45°  the  Nusselt  number  for  inclined  cylinders  is  greater  than  for  the 

same  array  in  cross-flow.  This  feature  is  almost  the  same  as  in  the  experimental  results  of 
Willins  and  Griskey  (1975),  who  measured  the  mass  transfer  from  a  single  cylinder  inclined 
relative  to  a  uniform  flow.  In  the  Willins  and  Griskey  experiment  the  group  Reo  Sc  was  of  the 
order  of  10^  (note  that  Reo  Sc  is  the  mass  transfer  equivalent  of  the  Peclet  number  in  forced 
convection  heat  transfer).  In  Fig.  1.6c,  the  Pr  =  7  curve  corresponds  to  Peo  =  70,  which  is  the 
highest  Peclet  number  illustrated  in  that  frame.  Moreno  and  Sparrow  (1987)  found  an  increase  in 
Nu  at  small  angles  P  in  flow  through  tube  banks  at  Reo  ~  10^  if  the  cylinders  were  in  line,  but 
not  if  they  were  staggered.  As  our  array  becomes  more  sparse,  it  more  closely  approximates  the 
in-line  configuration.  In  conclusion.  Fig.  1.6c,  Willins  and  Griskey's  and  Moreno  and  Sparrow’s 
results  suggest  that  the  increase  in  Nu(p)  above  Nu(0°)  in  a  small  p  range  is  a  high  Peclet  number 
effect  visible  in  in-line  arrays,  sparse  staggered  arrays  or  isolated  cylinders. 

Figure  1.7  shows  how  the  Nusselt  number  behaves  at  constant  Pr  and  constant  Rcd-  The 
curves  show  that  the  porosity  has  a  relatively  weak  effect  and  that  the  ratio  Nu(P)/Nu(0°)  is 
mainly  a  function  of  p.  Again,  the  exception  is  the  ^  =  0.95  curve  in  Fig.  1.7b,  which  shows  that 
the  porosity  plays  a  role  when  the  array  is  sufficiently  sparse. 

Figure  1.8  shows  that  the  Rej)  effect  is  such  that  Nu(p)/Nu(0°)  increases  when  Reo 
increases.  This  is  interesting  because  it  means  that  the  Nu(p)/Nu(0°)  curve  becomes  less  and  less 
like  Zukauskas's  curve  (Fig.  1.5)  as  Reo  increases.  This  trend  is  the  opposite  of  what  we  found 
in  the  results  for  the  pressure  drop,  where  the  AP(p)/AP(0°)  curve  approached  Zukauskas's  curve 
as  Rcd  increased  (Fig.  1.2a). 

Finally,  in  Figs.  1.9a,  b  we  show  that  when  the  Peclet  number  (Peo  =  Rcd  Pr)  is  held 
constant  (of  order  100),  the  ratio  Nu(p)/Nu(0°)  is  essentially  a  function  of  p  only.  Plotted  are 
two  curves  (Peo  =  70  and  100)  for  ([)  =  0.6  and  (j)  =  0.9,  and  a  single  point  (Rod  =  30,  Pr  =  3.3) 
for  <1)  =  0.6.  Comparing  this  finding  with  the  relative  pressure  curves  of  Fig.  1.2  leads  to  the 
conclusion  that  when  the  relative  pressure  curve  is  independent  of  Reo,  (namely,  when  Reo  ^ 
10),  the  ratio  Nu  (P)  /  Nu  (0°)  is  a  function  of  p  and  Peo,  and  not  of  Rej)  and  Pr  separately. 

The  dependence  of  Nu(P)  /  Nu(0°)  on  Pr  and  Reo  indicates  that  the  principle  of 
independence  does  not  hold  in  a  strict  sense.  It  does,  however,  provide  a  good  order-of- 
magnitude  estimate  for  the  effect  of  p.  That  the  principle  of  independence  does  not  hold  exactly 
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Figure  1.6  The  effect  of  angle  of  inclination,  Prandtl  number  and  porosity  on  the  Nusselt 
number  (Ren  =10). 


in  the  low  Rod  regime  is  reasonable,  since  this  principle  is  based  on  the  assumption  that  the  flow 
away  from  the  cylinder  surface  is  inviscid. 

The  most  important  aspect  of  the  heat  transfer  results  presented  in  this  section  is  that  at 
low  Reynolds  numbers  the  effect  of  cylinder  inclination  is  relatively  small,  regardless  of  the 
Peclet  number.  For  example,  if  the  p  effect  is  ignored  in  the  estimation  of  Nu  in  the  P  range  0'- 
60°,  the  error  in  Nu  is  less  than  20  percent. 

1.6  Experiments  for  Pressure  Drop 

We  verified  the  accuracy  of  our  numerical  results  by  measuring  in  the  laboratory  the 
pressure  drop  across  a  long  bundle  of  cylinders  with  several  inclinations  (P  =  0°,  45°,  60°)  and  at 
low  Reynolds  numbers  (8  <  Reo  <  50).  We  chose  to  test  the  pressure  drop  calculations  because 

they  represent  the  most  critical  part  of  the  results  described  until  now,  since  the  effect  of  cylinder 
inclination  is  much  greater  on  the  pressure  drop  than  on  the  heat  transfer  coefficient.  We  needed 
to  verify  that  our  computational  domain  captured  correctly  the  three-dimensional  flow  effects 
present  in  an  array  with  many  yawed  and  long  cylinders.  Earlier  experimental  studies  of  low- 
Reo  flows  through  cylinders  in  cross-flow  (Omohundro,  1949;  Bergelin  et  al.,  1949,  1950)  were 
limited  by  the  fact  that  the  arrays  were  too  short,  and  the  flows  were  not  fully  developed. 
Furthermore,  we  could  not  find  any  experimental  reports  on  low-Reo  flows  through  inclined 
cylinder  arrays. 

1.6.1  Apparatus.  The  main  features  of  the  experimental  apparatus  are  shown  in 
Fig.  1.10.  We  constructed  an  array  with  12  cylinders  across  and  115  cylinders  in  the  flow 
direction,  with  the  cylinders  arranged  in  the  configuration  indicated  as  orientation  1  in  Fig.  1.1. 
The  length  of  the  array  was  60  cm.  The  cylinders  were  made  of  0.07  inch  (1.8mm)  rubber  O- 
ring  stock.  They  were  7.4  cm  long  and  spaced  in  equilateral  triangles  with  side  length  of  6mm. 
When  the  cylinders  were  oriented  across  the  flow  (P  =  0°)  the  porosity  of  the  array  was  ([)  =  0.92. 
The  channel  height  to  cylinder  diameter  ratio  was  420:1.  In  order  to  change  the  fiber  angle,  the 
wall  of  the  lower  channel  was  shifted  longitudinally  and  upward,  such  that  the  total  fiber  length 
remained  the  same  (7.4  cm)  but  the  distance  between  the  top  and  bottom  of  the  channel  changed, 
and  the  fiber  angle  changed.  When  P  =  60°  the  channel  height  to  cylinder  ratio  was  210:1  and 
the  porosity  was  ({)  =  0.84. 

The  fiber  channel  was  connected  to  a  2  m  tall  tank  into  which  fluid  could  be  pumped, 
thereby  developing  a  pressure  head  to  drive  flow  through  the  array.  Six  manometers  were  spaced 
along  the  channel  so  that  pressure  measurements  could  be  made  in  the  region  of  fully  developed 
flow.  The  fluid  drained  out  of  an  opening  at  the  end  of  the  fiber  channel  into  a  settling  tank, 
from  which  it  was  siphoned  into  the  pumping  tank  and  re-pumped  into  the  head  tank. 
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Figure  1.10  The  main  components  of  the  experimental  apparatus. 
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The  pump  operated  at  constant  flow  rate  during  each  run,  so  that  the  height  of  the  fluid 
column  in  the  head  tank  adjusted  itself  until  the  flowrate  through  the  fiber  channel  matched  the 
flowrate  into  the  head  tank.  We  measured  the  flowrate  by  timing  the  filling  of  a  container  placed 
under  the  exit  from  the  fiber  channel.  The  fluid  was  a  solution  of  com  symp  and  water  (ratio  1:3 
by  volume)  with  a  kinematic  viscosity  between  7-8  x  10^  m^/s.  The  viscosity  depended  on  the 
temperature  and  the  concentration  of  com  syrap,  and  varied  from  day  to  day  but  was  stable 
during  a  series  of  consecutive  mns.  It  was  measured  after  each  mn  with  a  Cannon-Fenske  type 
viscometer. 

The  pump  was  a  variable  speed  Simer  shallow  well  pump.  The  manometers  were  glass 
tubes  with  4  mm  inside  diameter  and  Im  height.  The  head  tank  and  the  walls  of  the  fiber 
channel  were  made  of  Lexan.  The  fibers  (total  number  =  1380)  were  tied  on  the  back  side  of 
orifices  drilled  in  the  top  and  bottom  walls  of  the  fiber  channel,  and  sealed  with  RTV  silicone 
sealant.  The  fibers  were  pulled  tight  to  prevent  bending  during  operation.  No  bending  or 
vibration  of  the  fibers  was  observed;  however,  only  the  outermost  fibers  could  be  seen,  due  to  the 
density  of  the  fibers. 

1.6.2  Procedure.  The  pump  was  set  at  a  certain  flowrate.  After  the  pump  had  ran  for  a 
few  minutes,  to  let  the  bubbles  flow  out  and  to  allow  the  head  tank  to  reach  its  equilibrium 
height,  the  liquid  levels  in  the  manometers  were  read  and  recorded.  Then  the  fill  time  for  a  one- 
gallon  container  collecting  fluid  at  the  fiber  channel  exit  was  measured  five  times.  Finally,  the 
manometer  levels  were  read  again  to  make  sure  that  steady  state  had  been  achieved. 

Four  of  the  manometers  were  positioned  in  the  region  of  fully  developed  flow.  After  a 
few  mns  it  became  clear  that  it  was  going  to  be  unusual  to  produce  a  run  in  which  all  four 
manometers  were  bubble  free  at  the  time  of  measurement.  In  order  to  maintain  consistency 
between  runs,  the  reading  of  one  manometer  was  discarded  from  every  mn.  Pressure  was 
determined  by  the  relationship  P  =  pgh,  where  p  is  the  measured  fluid  density,  g  is  gravitational 
acceleration,  and  h  is  the  height  of  the  column  measured  in  the  manometer.  The  fluid  density 
was  measured  by  weighing  a  10  ml  fluid  sample.  A  95  percent  confidence  interval  was 
calculated  based  on  the  three  pressure  measurements.  Only  in  the  case  of  one  run  was  the 
precision  error  found  to  be  greater  than  5  percent,  and  the  mh  was  discarded.  All  other  mns  for 
which  all  of  the  necessary  measurements  could  be  made  were  included  in  the  reported  results. 

The  longest  and  shortest  measured  times  for  flowrate  were  discarded,  and  the  remaining 
three  times  were  averaged  to  find  the  flowrate  and  to  calculate  the  precision  error.  The  fluid 
velocity  derived  from  this  measurement  was  the  average  channel  velocity,  which  is  equivalent  to 
the  inlet  velocity  Uoo  used  in  the  numerical  simulations. 

1.6.3  Error  Analysis.  All  errors  were  calculated  using  2  standard  deviations  as  the  95 
percent  confidence  interval.  There  were  four  sources  of  precision  error  in  this  experiment 
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(random  errors  and  unsteadiness):  the  height  measurements  in  the  manometers  Ph,  the  kinematic 
viscosity  measurements  Py,  the  density  measurements  Pp,  and  the  flowrate  measurements  Pp. 
The  precision  errors  Ph  and  Pp  changed  from  run  to  run.  The  precision  error  for  density  (Pp/p  ~ 
0.07  percent)  was  negligible  compared  to  other  errors  and  did  not  vary  from  run  to  run. 
Repeated  tests  for  kinematic  viscosity  in  which  fluid  samples  were  taken  from  different  tanks  in 
the  system  indicated  that  Py/v  =  2.2  percent. 

The  bias,  or  the  fixed  errors  for  height  and  distance  measurements,  were  calculated 
assuming  that  our  ruler  was  accurate  to  within  ±  0.5mm.  Because  the  heights  being  measured 
changed  depending  on  flowrate,  Bh/h  varied  with  each  run,  from  0.5  to  4  percent.  The  fixed 
error  for  the  gallon  container  was  estimated  by  repeated  calibration  against  highly  accurate 
volumetric  measuring  containers,  and  was  found  to  be  1  percent.  Due  to  the  extremely  accurate 
containers,  scales  and  viscometer  calibration,  the  bias  errors  for  those  measurements  were 
negligible  compared  with  the  height  and  gallon  volume  errors.  The  non-negligible  bias  errors 
and  the  precision  errors  that  did  not  vary  from  run  to  run  are  listed  below: 

Py/v  Bx/x _ Bv/V _ Pp/p _ 

2.19%  0.5%  1%  0.07% 

The  95  percent  confidence  was  determined  by  combining  all  bias  and  precision  errors 
using  the  root- sum- square  method,  and  is  indicated  graphically  for  each  point  plotted  in  Figs. 
1.11  and  1.12.  The  data  scatter  falls  within  the  estimated  errors. 

1.6.4  Results.  Figure  1.11  shows  the  experimental  points  for  p  =  0°,  and  the 

corresponding  numerical  results  (the  solid  line).  Note  again  that  U  is  the  volume  averaged 
velocity,  equivalent  to  Uc®  of  the  numerical  simulations.  The  Reynolds  numbers  covered  by  the 
data  vary  from  10  to  50.  Our  numerical  solution  becomes  less  accurate  for  Rod  >  30.  A  possible 
explanation  for  this  is  that  the  experimental  system  enters  the  transition  to  the  turbulent  regime 
for  higher  Reo-  This  behavior  is  suppressed  in  our  numerical  simulation  by  the  assumption  of 
impermeable  side  walls  for  the  computational  domain.  These  boundary  conditions  do  not  allow 
the  flow  to  meander.  The  numerical  and  experimental  results  show  very  clearly  that  for  low  Reo 
the  flow  is  approximately  Darcy  flow,  and  that  the  numerical  simulation  predicts  the 
permeability  of  the  array  extremely  well. 

Figure  1.12a  shows  the  experimental  results  for  P  =  60°  compared  with  the  numerical 
solution  that  corresponds  to  the  experimental  geometry.  The  number  of  elements,  the  boundary 
conditions  and  general  configuration  used  in  the  simulation  were  the  same  as  those  used  in  our 
earlier  numerical  work  (sections  1.2  and  1.3).  An  adjustment  was  made  to  account  for  the 
specific  geometry  of  the  experimental  system.  Unlike  our  earlier  numerical  work,  in  the 
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Figure  1.12  Cylinders  inclined  at  60'  relative  to  the  cross-flow  position:  (a)  Comparison 
between  experimental  and  numerical  results,  and  (b)  Comparison  between  the 
porous  medium  model  (1.12)  and  the  experimental  results. 


experiment  the  porosity  is  a  function  of  the  angle  P,  and  the  fibers  are  no  longer  arranged  in 
equilateral  triangles  when  P  0.  The  two  curves  in  Fig.  1.12a  show  the  numerical  solution  for 
the  experimental  geometry  at  P  =  0°  and  60°.  The  agreement  between  the  experimental  and 
numerical  results  is  excellent. 

The  results  illustrated  in  Fig.  1.12a  cover  the  Ren  range  10-50.  In  the  P  =  60’  array, 
where  the  porosity  is  lower  than  when  P  =  0°,  the  Darcy  regime  covers  a  wider  Rep  range,  and 
our  numerical  results  agree  with  the  measurements  reasonably  well  over  the  entire  Rep  domain. 

Figure  1.12b  shows  the  experimental  data  for  P  =  60*  next  to  the  prediction  based  on  the 
porous  medium  flow  model,  equation  (1.12).  In  this  model  we  combined  the  permeability  from 
equation  (1.13)  with  the  angle  effect  documented  in  Figs.  1.2b  and  1.2d.  The  agreement  between 
the  proposed  model  and  the  experimental  data  is  excellent.  We  cannot  make  a  similar 
comparison  for  P  =  0°  because  it  falls  outside  the  Darcy  regime  (the  array  orientation  is  1,  and 
the  porosity  is  greater  than  0.9). 

Figure  1.13  shows  the  experimental  and  numerical  results  for  P  =  45°.  There  are  no  error 
bars  on  these  data  points  because  in  this  configuration  we  could  not  verify  the  precision  error; 
the  angle  was  such  that  the  entrance  and  end  region  effects  penetrated  to  two  of  the  central 
pressure  taps  that  were  supposed  to  be  in  the  fully  developed  flow  region.  Thus  only  two 
pressure  taps  were  in  the  fully  developed  region,  and  only  one  pressure  gradient  measurement 
was  possible.  The  pumping  mechanism  was  not  precise  enough  to  duplicate  exactly  the  flow 
conditions  for  two  successive  runs,  and,  therefore,  to  repeat  runs  while  holding  the  velocity 
constant  was  not  possible.  Instead  the  experiment  was  run  repeatedly  at  approximately  the  same 
flowrate.  If  repeated  pressure  gradient  measurements  indicated  that  the  first  run  was  not 
corrupted  by  an  unseen  bubble  or  other  anomaly  -  i.e.  if  two  runs  of  roughly  the  same  flowrate 
showed  roughly  the  same  pressure  drop  -  the  measurement  produced  by  the  first  run  was  plotted. 
The  consistency  of  the  data  plotted  in  Fig.  1.13  is  strong  evidence  that  the  measurements  are 
sufficiently  accurate,  and  they  agree  very  well  with  the  numerical  prediction. 

1.7.  Conclusions 

In  this  chapter  we  documented  the  pressure  drop  and  heat  transfer  characteristics  of  flows 
through  bundles  of  parallel  cylinders  when 

(a)  the  Reynolds  number  is  low, 

(b)  the  flow  is  hydrodynamically  and  thermally  fully  developed,  and 

(c)  the  cylinders  may  be  inclined  relative  to  the  flow  direction. 

These  three  aspects  had  not  been  documented  before,  yet  they  are  very  important  in  applications 
with  length  scales  smaller  than  the  scales  of  conventional  heat  exchanger  technology. 
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Figure  1.13  Comparison  between  the  experimental  and  numerical  results  for  the  array  inclined 
at  45'. 
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The  results  of  this  study  demonstrate  that  the  parametric  domain  represented  by  aspects 
(a)-(c)  requires  special  attention,  and  that  significant  errors  may  occur  if  the  available  large-scale 
(heat  exchanger-type)  information  is  extrapolated  to  smaller  scales.  For  this  reason,  a  useful 
direction  for  future  research  would  be  to  generate  more  data  for  pressure  drop  and  heat  transfer  in 
the  domain  (a)-(c).  Another  direction  would  be  to  widen  the  Reynolds  number  range  to  Reo 
values  as  high  as  100,  i.e.  to  cover  the  transitional  regime  where  the  flow  behind  each  cylinder 
meanders  (e.g.  Minakami  et  al.,  1993)  and  the  symmetry  assumed  in  Fig.  1  (top)  breaks  down. 

1.8  Notation 

A  wetted  area 

A  dimensionless  wetted  area 

^fiow  dimensionless  flow  cross-sectional  area 

Bh,v,x  bias  errors 

c  constant 

cp  fluid  specific  heat  at  constant  pressure 

D  cylinder  diameter 

F  dimensionless  mean  velocity,  equation  (1.18) 

H  height  of  computational  domain 

H  dimensionless  height 

K  permeability  of  porous  medium 

L  length  of  computational  domain 

L  dimensionless  length 

n  constant 

Nu  row  Nusselt  number 

p  pressure 

P  dimensionless  pressure 

Pr  Prandtl  number 

PhJF,v,p  precision  errors 
r  exponent 
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Rep 

S 


Reynolds  number 
distance  between  cylinder  centers 
T  temperature 

To  surface  temperature 

Too  inlet  fluid  temperature 

Ui  velocity  components 

u,v,w  velocity  components 

U,V,W  dimensionless  velocity  components 

Uoo  inlet  velocity 

W  width  of  computational  domain 

W  dimensionless  width 

x,y,z  cartesian  coordinates 

X,Y,Z  dimensionless  cartesian  coordinates 

a  fluid  thermal  diffusivity 

p  angle.  Fig.  1.1 

ATim  log-mean  temperature  difference 

9  dimensionless  temperature 

|i  viscosity 

v  kinematic  viscosity 

p  fluid  density 

(j)  porosity 
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Chapter  2 

FORCED  CONVECTION  FROM  A  SURFACE  COVERED  WITH 

FLEXIBLE  FIBERS 


2.1  Introduction 

In  this  chapter  we  document  the  fundamental  friction  and  heat  transfer  characteristics  of 
a  surface  covered  by  a  layer  of  flexible  fibers.  The  associated  forced  convection  phenomenon 
has  important  applications  in  several  fields,  for  example,  the  enhancement  of  heat  transfer, 
the  control  of  the  boundary  layer,  heat  and  mass  transfer  from  plant  canopies,  and  biological 
oceanagrophy,  where  suspension  feeders  and  plant  life  are  affected  by  the  local  flow  behavior. 
Our  own  interest  in  surfaces  covered  with  flexible  fibers  was  stimulated  by  questions  of  how 
to  maximize  the  insulation  efi'ect  of  hair  on  mammals. 

The  model  and  flow  configuration  described  in  this  paper  have  as  starting  point  the 
recent  work  on  modelling  convection  through  porous  media  in  contact  with  fluids.  In  natural 
convection,  the  interaction  between  the  flow  through  the  porous  medium  and  the  flow  of  the 
pure  fluid  was  documented  by  Poulikakos  (1986)  and  Sathe  et  al.  (1988).  Forced  convection 
was  documented  by  Vafai  and  Thiyagaraja  (1987)  and  Poulikakos  and  Kazmierczak  (1987). 
Relative  to  these  studies  and  our  own  work  on  surfaces  covered  with  hair,  the  present  study 
focuses  on  two  new  aspects: 

1 .  The  fibers  that  form  the  porous  layer  are  flexible,  and 

2.  The  external  flow  is  originally  perpendicular  to  the  surface. 

Aspect  (1)  deserves  scrutiny  because  when  the  fibers  are  flexible  the  interaction  between 
the  porous  medium  and  the  external  fluid  is  more  complex  than  when  the  fibers  are  rigid. 
The  flow  that  penetrates  the  porous  layer  can  change  the  local  properties  (directional  per¬ 
meabilities)  of  the  solid  matrix,  which  in  turn  influence  the  flow. 

Aspect  (2)  is  a  more  realistic  geometric  feature  of  the  flow  past  a  finite-size  body  covered 
with  fibers  (e.g.  the  body  of  a  mammal).  The  approaching  fluid  strikes  the  body  perpen¬ 
dicularly,  penetrates  the  fiber  cover  in  the  stagnation  region,  and  later  flows  parallel  to  the 
surface  and  around  the  body. 

2.2  Physical  Model 

The  outer  flow.  To  study  the  effects  of  fiber  cover  and  fiber  bending  we  selected  the 
convection  heat  transfer  configuration  shown  in  Fig.  2.1.  The  flow  enters  the  computational 
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domain  [y  =  yo)  with  uniform  velocity  {Uoo)  and  temperature  (Tc).  The  solid  impermeable 
wall  (y  =  0)  is  maintained  at  a  different  temperature  (Th).  The  flow  and  temperature  fields 
are  symmetric  about  x  =  0.  In  the  region  above  the  fibers  (y  >  H),  the  flow  is  governed 
by  the  Navier-Stokes  equations.  The  steady  state  conservation  of  mass  and  momentum  for 
incompressible  flow  are 


du  dv 
dx  ^  dy 

du  du  1  dp  ,  d^u  d^u , 

^dx^^dy  p  5x  5x2  5y2 


dv  dv  1  dp  .  d^v  d'^v . 
^5x'*'^5y  p  5y  ^  5x2  5y2 

with  the  following  boundary  conditions  in  the  plane  of  symmetry: 


(2.1) 

(2.2) 

(2.3) 


u  =  0,  =  0  at  X  =  0 

dx 

The  pressure  was  set  equal  to  an  arbitrary  constant  (p  =  0)  at  x  =  T  and  y  —  yo.  The 
velocity  components  («,  v)  were  matched  at  the  y  =  H{x)  interface  with  the  volume  averaged 
velocity  components  (u/,u/)  of  the  flow  through  the  layer  with  fibers. 


u  =  Uf,  v  =  Vf  at  y  =  H{x)  (2-5) 

The  selection  of  velocity  boundary  conditions  for  the  interface  between  a  fluid  and  a  porous 
medium  saturated  with  the  same  fluid  has  generated  considerable  amount  of  work,  which  is 
reviewed  in  Nield  and  Bejan  (1992).  According,  to  the  criteria  developed  by  Vafai  and  Tien 
(1981)we  do  not  need  to  include  the  Brinkman  term  in  our  porous  flow  model.  It  was  shown 
by  Vafai  and  Thiyagaraja  (1987),  however,  that  one  cannot  match  shear  stresses  across  the 
interface  without  the  inclusion  of  the  Brinkman  term.  If  one  does  not  match  shear  stresses 
then  one  has  to  account  for  the  possibility  of  a  slip  in  the  tangential  velocity  of  the  pure  fluid 
at  the  interface,  as  was  demonstrated  by  Beavers  and  Joseph  (1967).  In  this  problem  the 
flow  is  mainly  perpendicular  to  the  surface,  therefore  the  pressure  field,  which  drives  the  flow 
within  the  porous  layer,  is  not  a  strong  function  of  the  frictional  drag  at  the  surface.  The 
pressure  field  is  almost  completely  determined  by  the  decrease  in  the  vertical  momentum 
of  the  approaching  fluid.  The  result  of  this  is  that  the  problem  is  insensitive  to  the  degree 
of  slip  at  the  fluid-porous  interface.  Numerical  tests  indicated  that  if  the  no-slip  condition 
is  replaced  with  a  free-slip  condition,  the  resulting  change  is  never  greater  than  2  %  in  the 
frictional  drag  experienced  at  the  wall,  or  in  the  heat  transfer. 
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The  velocity  boundary  conditions  for  the  flow  out  of  the  computational  domain  (x  =  L) 
weved^ufdx'^  =  0  and  d^v / =  0.  The  insensitivity  of  numerical  results  to  the  specification 
of  outflow  boundary  conditions  has  been  demonstrated  by  numerous  researchers.  These 
conditions  were  compared  with  the  zero  stress  conditions  p  —  2nduldx  =  0  and  dvjdx  —  ^ 
at  X  =  L.  The  results  for  wall  friction  were  the  same  for  both  sets  of  boundary  conditions, 
but  the  zero  stress  condition  caused  velocity  oscillations  at  the  outlet. 

The  energy  conservation  equation  and  boundary  conditions  for  the  region  with  pure  fluid 
are 


dT  dT  fd^T  d'^T\ 


(2.6) 


T  =  To  at  y  =  yo  (2.7) 

T  =  Tf  at  y  =  H  -  (2.8) 

In  the  last  of  equations  (2.8),  the  fluid  temperature  is  matched  to  the  porous  medium  tem¬ 
perature  (Tj)  at  the  interface  between  the  two  regions.  Not  shown  is  the  energy  continuity 
equation,  in  which  the  heat  fluxes  are  matched  across  the  y  —  H{x)  interface. 


dT 

=  0  at  x  =  0, 
ox 


d'^T 

dx'^ 


=  0  at  X  =  L, 


The  layer  with  fibers.  It  was  assumed  that  the  flow  in  the  region  with  fibers  {0  <  y  <  H) 
is  in  the  Darcy  regime.  In  a  representative  elementary  volume  (r.e.v.)  (Nield  and  Bejan, 
1992)  within  this  region,  the  solid  matrix  is  a  bundle  of  parallel  equidistant  fibers  of  diameter 
D.  The  volume-averaged  flow  direction  almost  never  coincides  with  the  local  (r.e.v.)  fiber 
direction  because  the  outer  flow  is  not  uniform  and  the  fibers  can  bend.  When  the  flow 
through  the  fibers  is  perpendicular  to  the  fiber  direction,  the  permeability  may  be  estimated 
using 


D2 


nO.8 


(2.9) 


125  [(1  - 

In  the  opposite  extreme,  when  the  flow  is  parallel  to  the  fibers,  the  permeability  has  been 
correlated  by  Sparrow  and  Loeffler  (1959), 


^  _  -Ml  «^V2  .  . 

D2  -  16<^(1  -cl>)  ^  ^ 

This  correlation  is  valid  for  (j)  >  0-8  and  for  fibers  arranged  in  an  equilateral  triangular  array. 

By  analogy  with  the  irreversible  thermodynamics  of  heat  conduction  through  an  anisotropic 
medium  (Bejan,  1988),  the  Darcy  flow  equations  for  the  general  case  where  the  flow  is  not 
parallel  to  the  fibers  ^  0,  Fig.  2.2)  are 
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where 


Uj  = 

Kx  dp  Kxy  dp 
p  dx  p  dy 

(2.11) 

U/  = 

Rxy  dp  ^^y  dp 

p  dx  p  dy 

(2.12) 

A,  = 

Ks.  cos^  /?  +  A||  sin"*  ^ 

(2.13) 

II 

K±  sin^  P  +  A||  cos^  /3 

(2.14) 

II 

(A||  —  Ax)  sin  cos  /5 

(2.15) 

When  equations  (2.11,  2.12)  are  substituted  into  the  mass  conservation  equation  for  the 
fiber  region, 


duj  _  p 

dx  dy 

the  result  is  a  partial  differential  equation  for  the  pressure  field. 


(2.16) 


,dK, 

^  dx 


+ 


dKxy ^  ,dK^y 
dy  >dx^  ^dx-^^^  dx 


dy  dy  ^  dy"^  dxdy 


=  0 


(2.17) 


This  equation  allows  for  the  fact  that  the  permeabilities  vary  from  one  r.e.v.  to  another 
inside  the  fiber  region.  The  boundary  conditions  for  equation  (2.17)  are 


dp 

dx 

dp 

dy 


0 

0 


P  =  P{H) 


at 

X  =  0 

(2.18) 

at 

y  =  0 

(2.19) 

at 

x  =  L 

(2.20) 

and  the  condition  that  p  is  continuous  across  the  y  =  H  interface,  which  is  the  fourth 
boundary  condition  needed.  The  pressure  p{H)  is  provided  by  the  solution  to  the  pure  fluid 
part  of  the  flow  field.  The  constant  pressure  condition  (2.20)  is  the  better  of  two  conditions 
that  were  tried  at  x  =  X.  The  difficulty  associated  with  equation  (2.20)  is  that  it  forces 
Vf  =  Q  sX  X  =  L.  This  feature  causes  a  sharp  change  in  the  vertical  velocity  component 
across  the  y  =  H  interface  if  Rejr,  is  large.  The  alternative  to  equation  (2.20)  was  dpjdx 
=  constant  at  x  =  L.  This  was  superior  in  cases  where  the  fibers  did  not  bend,  but  when 
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the  fibers  bent  a  region  of  highly  negative  pressure  formed  around  x  =  L  and  y  =  0.  This 
was  not  acceptable  given  that  p  =  0  was  the  ambient  pressure,  so  the  constant-p  condition 
(2.20)  was  adopted.  It  must  be  said  that  either  condition,  constant  p  or  constant  dp/dx^ 
led  to  similar  results  for  friction  and  heat  transfer  over  the  wall  region  0  <  x  <  L:  The 
relative  error  between  the  two  sets  of  results  decreased  from  10%  at  T/if  =  2  to  less  than 
1%  at  T/ff  =  5.  In  these  Lf  H  ratios,  H  is  the  original  height  of  the  fiber  layer  before  it  is 
deformed  by  the  flow,  or  simply  H  evaluated  at  a;  =  0. 

We  now  turn  our  attention  to  the  thermal  aspects  of  the  model.  In  each  r.e.v.,  the 
temperature  is  represented  by  two  values,  the  temperature  of  the  solid  (the  fibers,  T,)  and 
the  temperature  of  the  interstitial  fluid  (T/).  The  conduction  of  heat  along  the  fiber  is 
described  by  the  unidirectional  fin  conduction  model 


+  =  o  (2.21) 

where  s  is  the  curvilinear  coordinate  measured  along  the  fiber  (s  =0  is  the  base).  In  each 
r.e.v.  cut  perpendicular  to  the  fiber  direction,  the  sum  of  all  the  fiber  wetted  perimeters  is 
Ps,  and  the  sum  of  all  the  fiber  cross-sections  is  A^.  For  each  local  angle  of  fiber  inclination 
/?,  the  ratio  Ps/Ag  is  a  constant  dictated  by  the  fiber  diameter  and  the  fiber  density  (or  the 
porosity  (f)  along  the  y  =  0  wall,  where  all  the  fibers  are  perpendicular  to  the  wall). 

The  equation  for  the  conservation  of  energy  in  the  interstitial  fluid  is 


UfdTf  VfdTf 
(f>  dx  ^  (f)  dy  ^ 


'd^Tf  d^Ti 


+ 


h  A 


ht 


pCp 


m  -  Tj) 


(2.22) 


5x2  Qy2 

where  Keu  is  the  r.e.v.  volume,  and  Aht! {(jA/rev)  ^^e  ratio  of  the  total  fiber-fluid  contact 
area  present  in  the  r.e.v.,  divided  by  the  volume  inhabited  by  fluid  in  the  r.e.v.  The  thermal 
diffusivity  a  refers  to  the  fluid  alone,  a  =  k({(f>Cp).  The  heat  transfer  coefficient  appearing  in 
equations  (2.21,  2.22)  was  estimated  based  on  the  low  Reynolds  number  correlation  developed 
based  on  the  work  described  in  Chapter  1,  where  it  was  found  that  in  the  range  0°  <  ^0  <  60° 
the  effect  of  /3  on  h  is  negligible. 


h  /U 

=  (8.46- 6.8.^)^  (2.23) 

This  correlation  is  valid  in  the  range  0.9  <  </>  <  0.95, 1  <  UooDja  <  30,  and  0.72  <  Pr  <  100. 
The  boundary  conditions  for  equations  (2.21,  2.22)  are 


Ts 

dTf 

dx 


Tj  =Tha.ty  -0 
0  at  X  =  0  and  x  —  L 
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(2.24) 

(2.25) 


Tf 

dZ 

ds 


=  T  dXy  =  H 

=  0  at  the  fiber  end 


(2.26) 

(2.27) 


Worth  noting  are:  the  outflow  condition  dTfjdx  =  0  a,i  x  =  L,  which  is  commonly  used  in 
numerical  studies  of  convection  in  porous  media  with  permeable  walls,  the  fluid  temperature 
continuity  across  the  y  =  H  interface,  equation  (2.26),  and  the  assumption  that  each  fiber 
is  slender  enough  so  that  the  heat  transfer  through  its  tip  can  be  neglected,  equation  (2.27). 

The  equations  and  boundary  conditions  described  until  now  were  nondimensionalized  by 
defining 


X  = 


ju.Uf) 

f/oo  ’ 


H 

uZl 


-  V  =  l-  P=  P 

H  pUl 


(2.28) 


(2.29) 


{T,Tf,Z)-Z 


(2.30) 

For  brevity,  we  omit  the  dimensionless  equations  and  boundary  conditions,  and  note  that 
they  contain  the  following  dimensionless  groups: 


Rei  = 


UocL 


Pr  = 


a 


-I 


hL 

Nu,=  -, 


Aht  — 


AhiL 


Bi  = 


hL 


Ps  = 


PsL 


(2.31) 


(2.32) 


The  fiber  shape.  The  third  component  of  the  model  contains  the  equations  needed  for 
calculating  the  shape  of  the  fibers.  This  component  is  based  on  the  observation  that  in  Darcy 
flow  the  pressure  gradient  in  the  fluid  is  balanced  by  the  forces  exerted  on  the  solid  matrix. 
We  assume  that  the  fibers  do  not  touch,  i.e.  that  each  fiber  acts  independently  of  all  other 
fibers.  Consider  a  small  volume  dxdyW  that  contains  at  least  one  r.e.v.,  where  W  »  L 
is  the  dimension  perpendicular  to  the  plane  of  Fig.  2.1.  The  components  of  the  total  force 
experienced  by  the  solid  parts  (fiber  segments)  found  in  this  volume  are 


F,  =  -^dxdyW,  Fy  =  -^dxdyW 

The  number  of  fibers  that  pass  through  this  volume  is 


n  =  (cos  fddx  +  sin  ^dy)W 


i-4> 


(2.33) 


(2.34) 
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where  (1  —  cf))  j {-k j i)  is  the  fiber  density  (number/area)  in  the  plane  perpendicular  to  the 
local  fiber  direction.  Dividing  equations  (2.33)  by  n  yields  the  components  of  the  local  load 
experienced  by  a  single  fiber.  These  forces  can  be  arranged  into  a  stress  vector  a:  the  normal 
and  tangent  forces  exerted  at  a  certain  point  (s)  on  the  fiber  are  obtained  by  performing  the 
scalar  product  between  a  and  the  unit  normal  and  unit  tangent  vectors. 

Each  fiber  was  divided  into  many  small  segments  {N),  typically  =  50.  The  forces  were 
modelled  as  acting  on  the  ends  of  each  segment.  Each  small  segment,  therefore,  behaved 
as  a  beam  undergoing  small  deflection.  The  normal  and  tangential  forces  were  transferred 
down  the  fiber  according  to  the  equations 

FNj  =  <7,  ■  iV,  +  {N,  ■  y+,)r;v,+i  +  (JV,  ■  fj„)FT,^,  (2.35) 

FTj  =  »i  ■  f,  +  (fj  ■  f,^,)FTi^,  +  (fj  ■  N,^,)FNi„  (2,36) 

where  the  unit  normal  N  is  defined  as  (cos — sin  y0j)  and  the  unit  tangent  T  is  defined 
as  (sin/?j,  cos/?j).  Following  the  method  of  Knight  and  Barret  (1980),  these  forces  can  be 
used  to  calculate  the  deflection  of  each  fiber  segment,  and  thereby  build  up  the  shape  of  the 
entire  fiber. 

Figure  2.2  illustrates  the  node  placement  relative  to  one  fiber  segment,  where  7j  is  the 
local  fiber  angle  at  node  j.  The  angle  ^  associated  with  the  jth  segment  is  the  angle  between 
the  y  axis  and  the  straight  line  segment  joining  the  j  and  j  +  1  nodes.  The  fiber  segments 

that  join  the  nodes  are  not  necessarily  straight.  The  mass  stations  that  would  be  present  in 

a  time- depen  dent  beam  analysis  are  not  necessary  because  in  this  study  we  are  concerned 
with  the  steady  state.  Only  the  elastic  elements  are  represented. 


A'fj  = 


f  —  Ajj  -b  7j-i 

(2.37) 

MjAs  FNj{Asy 

El  ^  2EI 

(2.38) 

Mj^i  +  AsFNj^i 

(2.39) 

In  equations  (2.38,  2.39)  As  is  the  length  of  the  fiber  segment,  E  is  the  elastic  modulus, 
I  =  7rD^/64  is  the  area  moment  of  inertia,  and  Mj  is  the  bending  moment  at  node  j. 
Equations  (2.35)-(2.39)  and  the  end  conditions  7  =  0  at  s  =  0,  and  =  0  at  the  free  end 
are  sufficient  for  calculating  all  the  'jjS.  To  calculate  the  position  of  each  node,  we  calculate 
the  end  deflection  of  each  segment. 
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(2.40) 


M,(As)2  FN^iAsf 

'  2EI  3EI 

The  position  vector  for  node  j, 


Pj  =  Asr(7,-,)  +  rfjiV(7,-.)  (2.41) 

where  T  and  N  are  the  tangent  and  normal  vectors  are  calculated  in  terms  of  the  angles  7j 
rather  than  13 j.  The  local  fiber  angle  can  be  calculated  for  each  segment, 


/3j-i  =  tan  ^ 


\Pix  ~  P{i-i)x / 


(2.42) 


where  Pjy  refers  to  the  y  position  of  the  jth  node,  and  Pj^  refers  to  the  x  position  of  the 
;th  node.  These  equations  were  also  nondimensionalized,  and  this  brought  to  light  a  new 
dimensionless  group  that  describes  the  relative  stiffness  of  the  fibers. 


S  = 


El 


(2.43) 


2.3  Numerical  Method 

The  equations  governing  the  flow  in  the  pure  fluid  region  were  solved  using  finite  dif¬ 
ferences  and  ADI.  The  pressure  field  was  solved  in  order  to  force  the  conservation  of  mass 
using  the  auxiliary  potential  method  (Fletcher,  1991).  The  grid  in  the  pure  fluid  region  was 
uniform  in  re,  but  variable  in  y,  so  that  a  large  number  of  grid  points  could  be  put  in  the 
boundary  layer  near  the  porous  interface.  Fig.  2.3.  A  shearing  transformation  was  used  to 
map  the  nonrectangular  physical  domain,  which  occurs  when  the  fibers  bend,  onto  a  rect¬ 
angular  computational  domain.  The  robustness  of  the  solver  was  increased  dramatically  by 
including  an  artificial  viscosity  that  was  non-zero  only  in  a  region  with  positive  v.  There 
is  no  artificial  viscosity  in  the  final  results  because  all  the  solutions  are  characterized  by 
negative  v  velocities. 

The  pressure  field  in  the  layer  with  fibers  was  also  solved  using  finite  differences  and 
ADI.  The  grid  in  the  porous  region  was  uniform  in  x  and  y  in  computational  space.  Fig.  2.3. 
Another  shearing  transformation  was  used  to  map  the  irregular  porous  domain  onto  a  rect¬ 
angular  computational  space. 

The  equations  for  the  fluid  temperature  were  solved  using  ADI  on  a  non-sheared  grid, 
with  the  same  grid  point  placement  as  the  fluid  equations  when  the  fibers  were  unbent.  The 
temperature  in  the  solid  was  determined  using  time  stepping  with  the  trapezoidal  rule  until 
convergence  was  reached. 
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x/L 


Figure  2.3:  Example  of  the  grid  used  {Rcl  =  1000,  Lj H  —  2,^  =  0.9,  H/D  —  20) 
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The  pure  fluid  velocity  solution  was  iterated  until  the  relative  error  (maximum  change 
in  velocity  divided  by  the  time  step)  was  less  than  0.1.  The  auxiliary  potential  solver  was 
then  iterated  until  its  relative  error  became  less  than  0.001.  Then  the  porous  pressure  field 
was  solved,  using  the  new  pressure  boundary  conditions  provided  by  the  auxiliary  potential 
solution.  The  porous  pressure  field  was  iterated  until  its  relative  error  was  less  than  0.0001. 
When  the  fibers  were  bent,  convergence  testing  indicated  that  a  relative  error  of  10“®  was 
necessary. 

The  cycle  continued  until  the  average  value  of  the  auxiliary  potential  function  was  less 
than  0.001  (which  indicates  that  continuity  is  being  satisfied).  All  of  these  tolerances  were 
determined  doing  convergence  testing.  In  all  cases,  dividing  any  of  the  above  convergence 
criteria  by  2  resulted  in  a  less  than  1%  change  in  the  calculated  friction  force  on  the  y  =  0 
surface.  Two  additional  accuracy  criteria  that  were  met  are  that  global  mass  conservation 
must  be  satisfied  to  better  than  1%,  and  that  mass  conservation  within  the  porous  layer  must 
be  satisfied  within  1%.  The  latter  criterion  was  relaxed  to  3%  when  the  fibers  underwent 
significant  deformation.  It  was  also  found  that  setting  j/o  =  4i  made  the  solution  insensitive 
to  further  increases  in  the  height  of  the  computational  domain. 

Similarly,  the  fluid  temperature  was  advanced  by  one  time  step  (the  size  of  which  was 
determined  by  stability  criteria)  and  then  the  solid  temperature  was  iterated  until  conver¬ 
gence  was  reached,  and  the  fluid  temperature  was  advanced  again.  The  system  generally 
conserved  energy  to  within  1%  if  the  relative  fluid  temperature  error  was  less  than  0.001.  In 
cases  where  the  total  conservation  of  energy  was  not  satisfied  to  within  3%,  error  bars  are 
attached  to  the  plotted  results. 

Grid  refinement  tests  were  performed.  The  final  grid  was  chosen  such  that  further  grid 
doubling  in  any  direction  resulted  in  a  less  than  1%  change  in  the  calculated  surface  friction. 
The  final  grid  for  Rcl  =  500  was  101  points  in  x  for  both  the  porous  and  non-porous  regions, 
81  grid  points  in  y  for  the  pure  fluid  region,  and  50  grid  points  in  y  for  the  porous  region.  An 
example  of  the  mesh  used  for  a  solution  in  which  the  fibers  are  bent  is  illustrated  in  Fig.  2.3: 
note  that  the  height  of  the  porous  layer  decreases  as  x  increases,  because  the  fibers  bend. 

2.4  Surface  with  Inflexible  Fibers 

Before  examining  the  effect  of  fiber  bending  on  friction  and  heat  transfer,  it  is  necessary 
to  understand  the  convection  mechanism  when  the  fibers  do  not  bend.  How  stiff  the  fibers 
(or  how  large  S)  must  be  in  this  limit  is  one  of  the  results  presented  in  the  next  section. 

The  independent  parameters  in  the  system  are  Rcl,  Pr,  the  porosity  of  the  fiber  layer  </>, 
the  aspect  ratio  T/if,  the  fiber  height  to  fiber  diameter  ratio  H/D,  and  the  ratio  of  solid  to 
fluid  conductivities  ks/k.  In  this  study  the  fluid  is  assumed  to  be  air,  Pr  =  0.72.  If  the  fibers 
are  hairs  in  air  then  kg/k  =  10.  We  considered  several  Rcl  values  between  500  and  2000, 
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several  porosities  between  0.9  and  0.95,  the  Hj D  ratios  20  and  50,  and  the  aspect  ratios 
LfH  =  1  through  5.  The  results  are  presented  in  terms  of  overall  quantities:  the  average 
skin  friction  coefficient. 


F;ix 

‘  \PUIL 


and  the  average  Nusselt  number 


where 


(2.44) 


“  where  4  =  (2,45) 

Figure  2.4  shows  that  the  average  skin  friction  coefficient  depends  mainly  on  Lj and 
is  roughly  inversely  proportional  to  Lj  H.  The  insensitivity  of  Cj  to  changes  in  Rei^  (p  and 
Hj D  is  particularly  evident  when  LjH  >  2,  although  minor  even  when  L/H  <2.  The  effect 
of  increasing  Hf  D  is  to  increase  Cf,  while  the  effect  of  increasing  (j)  is  to  decrease  Cf. 

The  results  for  heat  transfer  are  more  complicated.  Figure  2.5a  shows  the  effect  of  the 
porosity  and  the  Reynolds  number:  the  total  heat  transfer  rate  increases  with  both  4>  and 
Rei.  More  interesting  is  the  alternative  shown  in  Fig.  2.5b,  where  NulINul^i,  is  the  ratio 
between  the  actual  Nul  value  (Fig.  2.5a)  and  the  value  calculated  for  the  same  Rei  in  the 
limit  where  the  wall  surface  is  bare  (no  fibers,  or  LjH  oo).  Interesting  is  how  the  ratio 
Nui! Nui^b  compares  with  1.  Below  a  critical  porosity  (roughly  0.94  in  Fig.  2.5b),  the  fibers 
provide  an  insulation  effect,  and  Nui,  <  Nuifi-  Above  the  critical  porosity,  the  permeability 
of  the  porous  layer  is  sufficiently  high,  and  the  cold  flow  entering  from  above  is  sufficiently 
strong,  that  the  fibers  act  as  fins  surrounded  by  cold  fluid.  In  this  regime  the  fibers  augment 
the  heat  transfer  from  the  wall,  Nul  >  Nui^b.  These  two  extremes,  i.e.  the  fact  that  hair¬ 
like  fibers  provide  insulation  in  some  cases  and  augmentation  in  others,  confirm  the  main 
point  of  the  theoretical  work  or  surfaces  covered  with  hair  (Bejan,  1990). 

The  effect  of  Rcl  on  the  heat  transfer  is  illustrated  in  greater  detail  in  Fig.  2.6.  The 
heat  transfer  through  the  wall  with  fibers  increases  almost  proportionally  with  Rei.  This 
increase  is  steeper  than  when  the  wall  is  bare,  because  Nul! Nui^b  also  increases  with  L. 
In  conclusion,  the  faster  flow  (higher  Rei)  removes  more  heat  because  of  two  effects:  (1) 
cold  fluid  comes  closer  to  the  wall,  as  in  the  case  of  a  bare  wall,  and  (2)  the  fibers  provide  a 
finning  effect. 

Figure  2.7  shows  the  end  result  of  a  search  (Chapter  1)  for  a  way  to  correlate  the  numerical 
heat  transfer  results  for  surfaces  covered  with  stiff  fibers.  The  parameter  chosen  on  the 
abscissa  is  mLl H,  where  m  is  the  fraction  of  the  total  external  flow  (UooL)  that  penetrates 
into  the  fiber  layer. 


♦  Highest 
0  Lowest 


=  1.58(m4)°-^  -  0.051  (2.47) 

NuL,b  H'  ^  ^ 

approximates  the  data  with  a  mean  error  of  5.2  %. 

It  is  to  be  expected  that  equation  (2.47)  will  fail  in  the  limit  mL/ H  — >  0,  where  there  is 
no  flow  through  the  fiber  layer.  In  that  limit  the  Nusselt  number  reaches  a  minimum  value, 
NuL,Tnin,  which  Can  be  evaluated  based  on  the  parallel  thermal  resistance  model  (trapped 
fluid  in  parallel  with  fibers) 

+  -  (2'«) 

If  the  Nul  value  calculated  with  equation  (2.47)  is  less  than  NuL,mini  then  the  correct 
value  is  Nui^min  given  by  equation  (2.48).  To  illustrate  this  observation,  consider  the  left 
most  datum  (a  square)  plotted  in  Fig.  2.7:  that  point  corresponds  to  ^  =  0.9,  L/H  = 
2,i7/D  =  50  and  Rei,  —  500,  for  which  equation  (2.48)  yields  NuL^min  =  2.8.  This  means 
that  NuL,minlNuL,b  =  0.26,  which  is  greater  than  the  ratio  found  using  equation  (2.47), 
NuLfNuL^b  =  0.2.  The  actual  Nul  value  determined  numerically  is  2.9,  which  is  very  close 
to  NuL^rnin  =  2.8  given  by  equation  (2.48). 

To  calculate  the  heat  transfer  based  on  the  correlation  (2.47)  or  Fig.  2.7,  we  need  a  way 
to  predict  mLl H.  Figure  2.8  shows  the  correlation 

=  O.OmRe];'^  -  0.029  (2.49) 

C/,6 

where  C/,6  is  the  average  skin  friction  coefficient  for  the  bare  wall.  The  mean  error  between 
the  data  and  equation  (2.49)  is  18.6  %.  The  abscissa  parameter  is  the  Reynolds  number  used 
in  convection  through  porous  media  (Nield  and  Bejan,  1992),  namely  Re.p  = 

The  skin  friction  and  heat  transfer  correlations  for  the  bare-wall  limit  are  presented  in 
Fig.  2.9.  The  dashed  line  shows  an  empirical  correlation  developed  by  Sparrow  et  al.  (1979) 
from  experiments  measuring  mass  transfer  from  a  plate  two  and  a  half  times  as  wide  as  it 
was  long,  which  was  subjected  to  incident  flow  at  angles  from  25°  to  90°.  The  form  of  their 
correlation  is  Nul  =  O.dZORe^l^^ Pr^^^Rei,  where  L*  =  ^AfC,  A  is  the  plate  area,  and  C 
is  the  plate  circumference.  Sparrow  et  al.  found  that  by  using  Rcl*  in  the  correlation  they 
eliminated  almost  entirely  the  geometric  dependence  of  Nul-  Their  experiments  showed. 
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however,  a  small  geometric  dependence  for  the  leading  factor  (0.939  above):  this  factor 
increases  as  the  width  of  the  plate  increases.  In  their  experiments,  the  Nul  increase  from  a 
narrow  plate  (width/length  =  0.4)  to  the  wider  plate  was  9%.  This  change  is  consistent  with 
the  12%  difference  between  the  dashed  line  and  our  Nui,  results  (the  solid  line),  because  our 
results  are  for  a  plate  of  infinite  width. 

2.5  The  Effect  of  Fiber  Bending 

We  examined  the  effect  of  fiber  bending  by  varying  the  stiffness  parameter  S.  We  con¬ 
sidered  three  cases:  {Rcl  =  500,  LfH  —  l,(l>  —  0.9,  HjD  =  20),  {Rcl  =  1000,  L/H  =  2,<f)  = 
0.9, HjD  =  20)  and  {Rei  —  2000,X/if  =  ^,<f>  =  0.95, HjD  =  20).  The  selected  (f>  and  H/D 
values  are  such  that  a  significant  fraction  of  the  external  flow  penetrates  into  the  fiber  layer. 
The  characteristics  of  the  fiber  cover  on  the  LJ H  =  1  and  i/if  =  2  surfaces  are  identi¬ 
cal,  so  they  correspond  physically  to  two  surfaces  with  identical  covering,  but  with  different 
lengths,  i  and  i/2,  respectively.  The  maximum  fiber  bending  calculated  was  ^  =  20°  for 
i/if  =  1,  /?  =  40°  for  i/if  =  2  and  ^  =  38°  for  i/if  =  4. 

Figures  2.10  -  2.12  illustrate  the  effect  of  fiber  bending  on  the  skin  friction,  heat  transfer 
and  flow  penetration  into  the  fibers.  In  the  geometries  considered  the  transition  from  “no  S 
effect”  to  “substantial  S  effect”  is  fairly  abrupt.  More  flexible  fibers  (lower  S  values)  mean 
smaller  Cj,  Nui,  and  m  values. 

Figure  2.10  illustrates  the  effect  of  fiber  bending  on  skin  friction.  For  the  three  surfaces 
the  effect  is  significant.  The  decrease  in  Cj  from  the  case  of  no  bending  to  the  maximum 
bending  calculated  is  about  15%  in  all  cases  (specifically,  14%,  18%  and  15%  for  T/ff  =  1,2 
and  4,  respectively).  Figure  2.11  shows,  however,  that  the  effect  of  fiber  bending  on  Nul  is 
small.  The  total  decrease  in  Nul  is  5%,  8%  and  6%  for  the  three  surfaces.  From  the  result 
for  surfaces  with  inflexible  fibers  one  might  conclude  that  the  magnitude  of  the  penetrating 
flow  m  must  also  be  insensitive  to  the  fiber  bending.  Examination  of  Fig.  2.12,  however, 
shows  that  this  is  not  the  case. 

The  effect  of  fiber  bending  on  the  flow  penetrating  the  fiber  layer,  m,  is  shown  in  Fig.  2.12. 
The  effect  is  quite  large:  45%  for  L/H  =  2  and  43%  for  L/H  =  4.  According  to  the 
relationship  developed  for  predicting  Nui,  based  on  m  for  a  surface  covered  with  inflexible 
fibers,  equation  (2.47),  a  45%  decrease  in  m  should  result  in  about  a  21%  decrease  in  Nu^-, 
but  as  was  shown  in  Fig.  2.11,  Nui  actually  decreases  by  only  5-8%  while  m  decreases  by 
43  -  45%.  This  suggests  that  while  less  flow  enters  the  fiber  layer  as  a  whole,  the  flow  must 
be  penetrating  to  the  warm  surface  at  a  rate  roughly  independent  of  flber  bending. 

A  possible  explanation  for  this  is  that  when  the  fibers  bend  they  orient  themselves  in 
a  direction  closer  to  perpendicular  to  the  downward  coming  flow:  this  decreases  the  layer’s 
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Figure  2.11:  The  effect  of  fiber  stiffness  on  heat  transfer 


permeability.  At  the  same  time,  however,  the  spreading  of  the  fibers  causes  an  increase  in 
porosity,  which  causes  an  increase  in  permeability.  The  spreading  tends  to  occur  near  a;  =  0, 
and  the  maximum  fiber  bending  occurs  near  x  =  L  (Fig.  2.13).  This  results  in  a  deflection 
of  the  X  =  L  fluid  away  from  the  fiber  layer,  but  allows  fluid  into  the  small-x  region  at  the 
same  (or  even  a  higher)  rate  as  the  rate  that  occurs  when  the  fibers  do  not  bend. 

Figure  2.13  illustrates  the  progressive  bending  of  the  fibers  as  the  stiffness  S  decreases  for 
the  Lj H  =  4:  surface.  Each  “fiber”  line  drawn  on  the  fiber  surface  represents  many  actual 
fibers.  The  number  represented  per  line  depends  upon  and  H/ D.  It  is  clear  from  the 
figure  that  the  region  near  x  =  L  is  much  less  permeable  to  downward  flow  because  of  the 
fiber  angle  and  decreased  porosity.  The  minimum  porosity  for  the  Lj H  =  A  surface  was 
(j)  =  0.935,  which  corresponds  to  ^  =  0.95  when  the  fibers  are  straight. 

To  summarize.  Figs.  2.10  -  2.13  show  that  the  bending  of  the  fibers  begins  to  have  an 
effect  when  the  stiffness  number  S  drops  below  a  certain,  critical  level.  We  recorded  this 
effect  quantitatively,  by  first  defining  the  critical  stiffness  number  Sc  as  the  S  value  where 
m  drops  to  90%  of  its  value  for  inflexible  fibers. 


m[S  =  Sc)  =  0.9  m(5'  — >  oo)  (2.50) 

We  then  determined  numerically  the  m{S)  dependence  for  a  series  of  12  cases,  of  which  only 
three  are  illustrated  in  Fig.  2.10.  The  resulting  Sc  values  can  be  correlated  as 


5'.= 


C 

\-4> 


(2.51) 


for  which  the  C  coefficient  is  reported  in  Fig.  2.14.  The  analytical  form  of  equation  (2.51)  was 
derived  based  on  equation  (2.38)  and  the  hypothesis  that  the  critical  stiffness  Sc  corresponds 
to  a  critical  fiber  angle  7c.  Figure  2.14  shows  that  in  the  parametric  domain  considered  the 
C  coefficient  does  not  depend  systematically  on  Rei.  The  most  we  can  say  is  that  C  is 
a  number  of  order  0.4.  That  a  correlation  of  type  (2.51)  should  hold  is  suggested  also  by 
Fig.  2.4,  which  implies  that  the  dimensionless  pressure  field  is  largely  insensitive  to  Rei,. 

Setting  C  =  0.4  in  equation  (2.51)  allows  us  to  predict  Sc  to  within  a  factor  of  2  for  all  the 
cases  tested.  We  expect  this  relationship  to  hold  for  high  Re^  flows  where  the  viscous  effects 
on  the  pressure  field  will  continue  to  be  negligible.  Equation  (2.51)  should  break  down  at 
sufficienty  low  Rei,  values,  because  in  that  direction  viscous  effects  become  important.  Low 
Rei  flows,  however,  are  also  characterized  by  little  or  no  penetrating  flow  in  the  fiber  layer, 
which  makes  these  flows  of  considerably  less  interest. 
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2.6  Concluding  Remarks 

In  this  chapter  we  investigated  the  interaction  between  a  solid  surface  and  an  external 
flow  when  the  solid  is  covered  with  a  layer  of  fibers.  The  objective  was  to  document  the  way 
in  which  the  fiber  layer  properties  affect  the  overall  heat  transfer  and  friction  characteristics 
of  the  surface. 

In  the  first  part  of  the  study  we  considered  the  limit  in  which  the  fibers  are  inflexible. 
The  results  presented  in  Figs.  2. 4-2. 9  and  in  the  text  show  that  the  fiber  layer  can  alter 
significantly  the  wall  characteristics,  especially  the  overall  heat  transfer.  A  peculiar  effect 
is  due  to  the  porosity:  at  porosities  lower  than  a  critical  value,  the  fiber  cover  acts  as  an 
insulation,  while  at  higher  porosities  the  fibers  augment  the  heat  transfer  (Fig.  2.5b).  We 
showed  that  the  effect  of  the  fiber  cover  can  be  correlated  in  terms  of  the  flow  fraction  that 
penetrates  into  the  fibers  (Fig.  2.7,  2.8). 

In  the  second  part  we  examined  the  bending  of  the  fibers,  as  they  interact  with  the 
flow  that  penetrates  the  fiber  layer.  The  fiber  bending  is  governed  by  the  nondimensional 
stiffness  number  S.  We  found  that  when  S  decreases  below  a  critical  value  Sc,  there  is  a 
sudden  drop  in  wall  friction,  heat  transfer  and  penetrating  flow  fraction  (Fig.2.10-2.12).  The 
least  sensitive  to  changes  in  S  is  the  overall  heat  transfer  rate  (Fig.  2.11).  We  showed  that 
the  critical  stiffness  number  Sc  values  determined  numerically  can  be  correlated  cf.  equation 
(2.51)  and  Fig.  2.14. 

2.7  Notation 

A  aspect  ratio,  Lj H 
Aht  heat  transfer  area  in  one  r.e.v.,  m? 

Aht  dimensionless  heat  transfer  area 
As  sum  of  fiber  cross-section  in  one  r.e.v., 

Bi  Biot  number 

Cp  specific  heat  at  constant  pressure,  kJf(kgK) 

C  coefficient  for  the  critical  stiffness  number 
correlation  (2.51),  Fig.  2.14 
C f  average  skin  friction  coefficient 
dj  deflection  of  fiber  segment,  m 
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fiber  diameter,  m 

modulus  of  elasticity,  N j m? 

force  components,  N 

heat  transfer  coefficient 

height  of  unbent  fibers,  m 

area  moment  of  inertia, 

fluid  thermal  conductivity,  Wl{m  ■  K) 

fiber  thermal  conductivity,  Wl{m  ■  K) 

permeabilities,  w? 

half-length  of  plane  wall,  m 

penetrating  flow  fraction 

moment,  N  ■  m 

number  of  fibers  in  one  r.e.v. 

number  of  fiber  segments 

average  Nusselt  number 

pressure,  Njm? 

sum  of  wetted  perimeters  in  one  r.e.v.,  m 

dimensionless  sum  of  perimeters  in  one  r.e.v. 

dimensionless  pressure 

position  vector,  m 

Prandtl  number,  v  J  a 

heat  flux,  Wim? 

representative  elementary  volume 
Reynolds  number 
porous  medium  Reynolds  number 
curvilinear  coordinate,  m 
stiffness  number 
critical  stiffness  number 
temperature,  K 
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u,  V  velocity  components,  m/s 
U,  V  dimensionless  velocity  components 
Uoo  approach  velocity,  m/s 
Vrev  volume  of  one  r.e.v., 

W  width  of  plane  wall,  m 
x,y  coordinates,  m 
X,  Y  dimensionless  coordinates 
a  fluid  thermal  diffusivity, 

^  fiber  angle.  Fig.  2.2 
7c  critical  fiber  angle 
7j  local  fiber  angle  at  node  j 
9  dimensionless  temperature 
fi  viscosity,  kgl{s  •  m) 

V  kinematic  viscosity,  m^/s 
p  fluid  density,  kglm^ 
a  stress  vector,  Njinn? 

(f>  porosity 
(  )b  bare  wall 

(  )c  cold  fluid 

(  )h.  hot  wall 

(  )  f  fluid  in  the  fiber  layer 

(  )s  solid  matrix  (fibers) 

(  )”  per  unit  area 

(  ~  )  dimensionless 
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Chapter  3 


THE  OPTIMAL  SPACING  BETWEEN  HORIZONTAL  CYLINDERS  IN 
A  FIXED  VOLUME  COOLED  BY  NATURAL  CONVECTION 

3.1  Introduction 

In  this  chapter  we  report  the  results  of  a  theoretical,  numerical  and  experimental  study 
of  how  to  determine  the  optimal  spacing  for  horizontal  cylinders  in  an  array  with  natural 
convection  heat  transfer.  The  volume  occupied  by  the  array  is  fixed.  The  number  of  cylinders 
in  the  array,  or  the  spacing  between  cylinders  of  fixed  diameter,  can  vary.  The  optimal 
spacing  reported  in  this  paper  corresponds  to  the  maximum  overall  heat  transfer  (or  thermal 
conductance)  between  the  array  and  the  surrounding  fluid. 

The  optimal  spacing  question  is  both  important  and  timely.  It  is  important  because  of 
its  obvious  implications  in  the  design  of  heat  exchangers,  surfaces  with  horizontal  pin  fins, 
and,  generally,  the  cooling  by  natural  convection  of  a  space  of  fixed  size  (e.g.  electronic 
package).  The  question  is  timely  in  view  of  the  large  volume  of  research  that  has  been 
devoted  to  arrays  of  horizontal  cylinders  with  natural  convection  on  the  outside.  This  body 
of  work  was  reviewed  on  several  occasions  (e.g.  Guceri  and  Farouk,  1985;  Sadeghipour  and 
Asheghi,  1994)  and  will  not  be  reviewed  here.  Most  of  this  work  delt  primarily  with  the 
detailed  interaction  between  adjacent  cylinders,  not  with  the  overall  performance  (thermal 
conductance)  of  the  entire  array. 

The  existence  of  an  optimal  spacing  for  maximum  heat  transfer  was  noted  by  Sparrow 
and  Vemuri  (1985)  in  an  experimental  study  of  an  array  with  a  large  number  of  horizontal  pin 
fins.  The  heat  transfer  was  by  combined  natural  convection  and  radiation.  The  maximum 
exhibited  by  the  overall  heat  transfer  rate  was  very  shallow  and  corresponded  to  using  an 
array  with  approximately  35  pin  fins. 

Related  aspects  of  this  subject  were  considered  by  Tokura  et  al.  (1983)  and  Sadeghipour 
and  Asheghi  (1994).  Tokura  et  al.  reported  an  optimal  spacing  for  a  vertical  column  of 
horizontal  cylinders  confined  by  two  vertical  plates.  It  is  not  at  all  clear,  however,  that  their 
result  corresponds  to  maximum  heat  transfer.  Their  recommendation  to  use  6  cylinders  in 
the  array  appears  to  be  a  trade-off  between  transferring  more  heat  and  using  more  hardware 
(cylinders). 

Sadeghipour  and  Asheghi  reconsidered  the  optimal  spacing  question  in  Tokura  et  al.  ’s 
configuration.  Their  experimental  results  suggest  a  shallow  maximum  in  the  variation  of 
the  overall  heat  transfer  rate  with  the  cylinder-to-cylinder  spacing.  The  final  correlation  for 
this  relationship,  however,  does  not  indicate  a  maximum  with  respect  to  spacing:  the  overall 
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heat  transfer  increases  monotonically  with  the  spacing,  and  reaches  its  highest  values  (or 
plateau)  when  the  Sj D  ratio  exceeds  approximately  15. 

The  absence  of  a  correlation  for  the  reported  optimal  spacings,  and  the  inability  of  pre¬ 
dicting  optimal  spacings  in  new  configurations  defined  the  objectives  of  the  present  study. 
To  begin  with,  we  focused  exclusively  on  heat  transfer  by  natural  convection  (without  radi¬ 
ation).  Our  objective  was  not  only  to  determine  the  optimal  spacings  for  various  situations, 
but  also  to  correlate  our  results  in  a  compact  dimensionless  relation  valid  for  any  array  with 
laminar  natural  convection. 

To  achieve  this  objective  we  conducted  the  study  in  three  phases.  In  the  first,  we  devel¬ 
oped  a  pure  theory  to  prove  the  existence  of  the  optimal  spacing,  and  to  reveal  the  proper 
dimensionless  groups  and  analytical  form  of  the  optimal-spacing  correlation.  In  the  second 
phase,  we  simulated  the  natural  convection  flow  and  temperature  fields,  and  varied  the  spac¬ 
ing  (or  number  of  cylinders)  of  the  array  to  generate  optimal-spacing  data  to  be  correlated 
based  on  the  theory.  Finally,  in  the  third  phase  of  the  study  we  experimented  with  several 
arrays  of  cylinders  that  occupied  the  same  volume.  The  purpose  of  these  experiments  was  to 
determine  the  optimal  spacing  through  direct  heat  transfer  and  temperature  measurements, 
and  to  test  in  this  way  the  accuracy  of  the  numerical  phase  of  the  study. 

3.2  Theoretical  Results 

We  begin  with  a  theoretical  argument  to  demonstrate  that  there  exists  an  optimal 
cylinder-to-cylinder  spacing  for  maximum  heat  transfer,  and  to  identify  the  proper  dimen¬ 
sionless  groups  needed  to  correlate  the  optimal  spacing  results  determined  more  accurately 
(numerically  and  experimentally). 

Consider  the  bundle  of  horizontal  cylinders  shown  in  Fig.  3.1.  The  overall  dimensions 
of  the  bundle  [H,  L,  W)  and  the  cylinder  diameter  [D)  are  fixed.  Natural  convection  heat 
transfer  (g)  occurs  between  the  cylinder  surfaces  (^^u)  and  the  surrounding  fluid  reservoir 
{Too)-  In  the  following  analysis  we  address  the  question  of  how  to  select  the  number  of 
cylinders  in  the  bundle,  or  the  cylinder-to-cylinder  spacing  (5),  such  that  the  overall  thermal 
conductance  between  the  bundle  and  the  ambient,  qKT^  —  Too)  is  maximized.  For  the 
sake  of  conciseness  we  assume  that  the  cylinders  are  staggered,  and  that  their  centers  form 
equilateral  triangles.  Other  array  types  can  be  treated  similarly.  The  problem  is  similar 
to  that  of  determining  the  optimal  number  of  plates  in  a  stack  cooled  by  natural  or  forced 
convection. 

When  the  spacing  S  and  the  Rayleigh  number  are  suflaciently  large  each  horizontal 
cylinder  is  coated  by  a  distinct  boundary  layer,  and  the  surrounding  fluid  is  at  the  tem¬ 
perature  Too-  We  are  assuming  that  {H,W)  »  {D  -f  S),  and  that  Rao  »  1,  where 
Raf)  —  g^D^{Tyj  —  Too)l{oiv).  The  heat  transfer  from  one  cylinder  is 
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qi  =  ^NuD7rDL{T^  -  T^) 

where  the  overall  Nusselt  number  is  (cf.  Morgan,  1975,  10“*  <  Ran  <  10’^) 


(3.1) 


Nud  =  0,48fiai(* 

The  total  number  of  cylinders  in  the  bank  of  cross-sectional  area  77  x  IT  is 


(3.2) 


_  HW 

[S  -b  DY  cos  30° 

therefore  the  total  heat  transfer  from  the  bank  \s  q  =  nqi  or 


(3.3) 


q,„,es  =  1.74  (3.4) 

This  result  shows  that  when  the  spacing  is  large,  the  overall  thermal  conductance  qj [Tyj—Too] 
decreases  as  S  increases. 

Consider  now  the  opposite  extreme  when  the  cylinders  almost  touch,  and  the  flow  is 
almost  cut  off.  In  this  limit  the  temperature  of  the  coolant  that  exits  slowly  through  the 
upper  plane  of  the  bundle  {L  x  IT)  is  essentially  the  same  as  the  cylinder  temperature  T^. 
The  heat  transfer  from  the  bundle  to  the  coolant  is  equal  to  the  enthalpy  gained  by  the 
coolant,  q  =  mc.p{Tw  —  Too),  where  m  is  the  mass  flow  rate  through  the  L  x  W  plane. 

To  obtain  an  order-of-magnitude  estimate  for  the  flow  rate,  we  note  that  m  is  composed 
of  several  streams  [total  number  =  W/{S  +  D)],  each  with  a  cross-sectional  area  5  x  T  in  the 
plane  of  one  horizontal  row  of  cylinder  axes.  The  thickness  of  the  channel  traveled  upward 
by  each  stream  varies  between  a  minimum  value  (5)  at  the  row  level,  and  a  maximum  value 
at  a  certain  level  between  two  rows.  The  volume-averaged  thickness  of  one  channel  is 


5  =  5  +  T»  -  0.907-—  (3.5) 

O  T" 

however,  we  may  adjust  this  estimate  by  using  1  in  place  of  the  factor  0.907  to  account  for 
the  fact  that  the  channel  closes  (i.e.  the  flow  must  stop)  when  the  cylinders  touch  (5  =  0): 


S  =  S 


S  +  2D 
S  +  D 


(3.6) 


When  S  is  sufficiently  small,  the  flow  rate  through  each  channel  of  cross  sectional  area 
SL  and  flow  length  H  is  proportional  to  the  pressure  difference  that  drives  the  flow.  The 
pressure  difference  is  AT  =  pgHY{T.^  -  Too),  or  the  difference  between  the  hydrostatic 
pressures  under  two  T-tall  columns  of  coolant,  one  filled  with  Too  fluid,  and  the  other  with 
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(a)  (b) 


Figure  3.1:  Bundle  of  horizontal  cylinders  immersed  in  a  quiescent  fluid  (a),  and  detail  of 
one  of  the  channels  traveled  by  the  fluid  (b). 
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Tyj  fluid.  The  mean  velocity  of  the  channel  flow,  U,  can  be  approximated  using  the  Hagen- 
Poiseuille  solution  for  flow  between  two  parallel  plates  (spacing  S.,  flow  length  i7), 


S^AP 

l2nH 


(3.7) 


The  total  flow  rate  through  the  bundle,  m{pUSL)  ■Wl{SpD),  leads  to  the  total  heat 
transfer  rate  mCp{T.w  —  Too),  which  can  be  summarized  as 


^smalls  —  22fl)3( p  £)'^  ^{Tw  Too)R(^D  (^■^) 

The  key  feature  of  this  estimate  is  that  when  S  is  small  and  decreases,  the  thermal  conduc¬ 
tance  ql{Tyj  —  Too)  decreases  as  S^. 

Figure  3.2  summarizes  the  trends  uncovered  so  far.  The  actual  thermal  conductance 
would  be  represented  by  the  solid  curve  sketched  in  the  figure.  The  peak  of  this  curve 
corresponds  to  a  spacing  (Sopt)  that  can  be  approximated  by  intersecting  the  two  asymptotes, 
qiargeS  =  ^smaiis-  The  result  of  eliminating  q  between  equations  (3.4)  and  (3.8)  is 


Sopt  2  +  Sopt!  D 


^n1/3 


This  relation  is  plotted  in  Fig.  3.3,  which  shows  that  Sopt! D  is  almost  proportional  to  the 
group  {Hj In  other  words,  a  simpler  alternative  to  the  order  of  magnitude 
estimate  obtained  in  equation  (3.9)  is 


or 


RapI* 


(3.10) 


(3.n) 

where  Ran  =  —  Too)l{ocv).  Equation  (3.11)  shows  that  the  optimal  spacing  is 

approximately  proportional  to  which  means  that  it  is  almost  insensitive  to 

changes  in  the  cylinder  diameter. 


3.3  Numerical  Results 


In  this  section  and  the  next  we  report  numerical  and  experimental  results  for  the  optimal 
cylinder-to-cylinder  spacing.  The  objective  of  these  empirical  phases  of  our  study  was  to  test 
and  improve  the  accuracy  of  the  correlation  determined  theoretically,  equations  (3.9)-(3.11). 
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The  natural  convection  and  heat  transfer  in  an  array  of  horizontal  cylinders  was  simulated 
numerically  by  focusing  on  a  vertical  channel  formed  between  two  adjacent  rows  of  cylinders, 
Fig.  3.4.  It  was  assumed  that  the  flow  is  laminar,  such  that  there  is  no  exchange  of  fluid 
and  energy  between  adjacent  channels.  The  regime  is  laminar  when  Rajj  <  lO^Pr  (Bejan 
and  Lage,  1990),  because  H  is  the  relevant  vertical  dimension  of  the  channel  with  boundary 
layer  flow. 

The  mass,  momentum  and  energy  equations  were  simplified  in  accordance  with  the  as¬ 
sumptions  of  two-dimensional  steady  state,  nearly  constant  properties,  and  Boussinesq  ap¬ 
proximation  in  the  buoyancy  term  of  the  momentum  equation  for  the  vertical  (y)  direction: 


u 


u 


+ 


+ 


du 
dx 
du 
dx 
dv 
dx 
dT 

+ 


1^  =  0 
dy 

du  1  dp  ^2 

v—  = - ^  +  uV^u 

ay  p  ox 


dv 


1  dp 


v^  =  aV^r 

dy 


(3.12) 

(3.13) 

(3.14) 

(3.15) 


where  =  d'^ jdx'^  -f  The  horizontal  and  vertical  velocity  components  are  u  and 

V.  The  origin  of  the  Cartesian  frame  (x,?/)  is  located  in  the  bottom  left  corner  of  the 
computational  domain.  The  computational  domain  contains  the  actual  channel  of  height 
HjD  (fixed),  plus  an  inlet  (bottom)  section  and  an  outlet  (top)  section.  Accuracy  tests 
showed  that  when  the  inlet  length  is  2.5  D  and  the  outlet  length  3Z),  the  calculated  heat 
transfer  from  the  channel  is  insensitive  (with  changes  less  than  1  percent)  to  further  doubling 
of  the  inlet  and  outlet  lengths. 

The  flow  boundary  conditions  were:  zero  normal  stress  and  vertical  flow  (u  =  0)  at  the 
inlet  to  the  computational  domain  (x  =  0);  free  slip  and  no  penetration  at  the  fluid  interfaces 
(planes  of  symmetry)  between  two  consecutive  cylinders;  no  slip  and  no  penetration  at  the 
cylinder  surfaces,  and  free  slip  and  no  penetration  on  the  vertical  boundaries  of  the  inlet 
section.  It  is  important  to  note  that  the  forcing  of  free  slip  and  no  penetration  conditions 
on  the  vertical  boundaries  of  the  outlet  section  would  induce  an  artificial  acceleration  of 
the  fluid  (updraft,  chimney  effect)  through  the  cylinder-to-cylinder  channel.  To  avoid  this 
effect,  a  zero  stress  inlet  condition  was  specified  along  one  of  the  sides  of  the  outlet  section, 
specifically,  on  the  side  opposite  the  topmost  cylinder  (Fig.  3.4). 

The  temperature  boundary  conditions  were  T  =  Tu,  on  the  cylinder  surfaces,  and  T  =  T^o 
at  the  bottom  end  of  the  computational  domain.  The  remaining  portions  of  the  boundary 
of  the  computational  domain  were  modelled  as  adiabatic. 
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Figure  3.2:  The  optimal  cylinder-to-cylinder  spacing  for  maximum  thermal  conductance,  as 
the  intersection  of  the  large-S  and  small-S  asymptotes. 
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Equations  (3.12)  -  (3.15)  were  nondimenionalized  by  defining  the  dimensionless  variables 


{x,y) 
D  ’ 

r-T. 


([/,!/)  = 


(u,u) 


p 


{aj  D){RaDPryi^ 

pD‘^ 


(3.16) 


(3.17) 


Tw  —  Too'  ^  fia{RaDPry^‘^ 

Numerical  solutions  were  generated  for  Pr  =  0.72.  The  system  of  equations  was  solved 
on  a  Cray-YMP  using  a  finite  element  package  (FIDAP,  1991).  Grid  refinement  tests  showed 
that  40  nodes  per  D  in  x  and  y  were  necessary  for  Ra^  =  10^,  such  that  the  further  doubling 
of  the  number  of  nodes  caused  a  change  of  less  than  2  percent  in  the  heat  flux  from  each 
cylinder. 

The  accuracy  of  the  numerical  method  was  checked  further  by  placing  a  single  half 
cylinder  in  the  computational  domain,  and  calculating  the  overall  Nusselt  number  Nud- 
The  values  obtained  were  Nud  =  3.20  for  Rao  =  10^,  and  Nud  =  5.10  for  Rajj  —  10“*. 
These  values  agree  within  6  percent  with  the  large  body  of  empirical  results  compiled  by 
Morgan  (1975),  and  fall  between  the  error  bars  indicated  by  Morgan  for  the  correlation  listed 
in  equation  (3.2). 

For  each  array  geometry  {H/D^  S/D)  the  solutions  were  generated  by  solving  for  Rud  = 
10^  first,  and  then  using  that  solution  as  initial  condition  for  the  next  higher  Rao.  The 
computational  time  depended  strongly  on  the  array  geometry  and  Rayleigh  number. 

Figure  3.4  shows  a  sample  of  the  flow  and  temperature  fields  calculated  for  if/D  =  10 
and  Raf)  =  10"*.  The  streamline  and  isotherm  patterns  are  arranged  in  a  way  that  illustrates 
our  search  for  the  optimal  cylinder-to-cylinder  spacing  of  an  array  in  a  fixed  volume.  We 
fixed  the  height  of  the  array  (W),  and  decreased  the  spacing  {S)\  we  did  this  discretely  (in 
steps),  each  time  by  adding  one  more  to  the  number  of  cylinders  {nfi)  in  the  channel  of  height 
H. 

During  this  sequence  we  monitored  the  total  heat  transfer  from  all  the  half- cylinders,  to 
see  how  it  responds  to  changes  in  S.  A  sample  is  shown  in  Fig.  3.5,  where  the  ordinate 
parameter  represents  the  dimensionless  heat  transfer  volumetric  density. 


HLWk{T^-T^) 

where  q  is  the  total  heat  transfer  from  an  array  of  fixed  volume  HLW.  Figure  3.5  shows  that 
there  is  an  optimal  spacing  for  maximum  heat  transfer.  The  Sopt/D  values  were  determined 
by  fitting  the  highest  three  q  points  with  a  parabola,  and  solving  dq/d{S/D)  =  0.  The 
best  practical  spacing,  of  course,  corresponds  to  the  number  of  cylinders  n  (an  integer)  that 
maximizes  the  total  heat  transfer  rate  from  the  given  space. 
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Figure  3.4:  The  computational  domain,  and  the  effect  of  cylinder-to-cylinder  spacing  on  the 
flow  and  temperature  fields  {H/D  =  10,  Rao  —  10^,  Pr  =  0.72).  The  streamliners  are  on 
the  left,  and  isotherms  on  the  right. 


The  resulting  SoptjD  data  are  reported  in  Fig.  3.3.  They  are  correlated  very  nicely  in 
the  manner  anticipated  in  section  3.2,  namely  by  using  the  group  {Hj Dyl^Ra  on  the 

abscissa.  The  optimal  spacings,  however,  are  consistently  2.5  times  larger  than  the  values 
calculated  based  on  equation  (3.9).  The  function  of  type  (3.10)  that  fits  the  numerical  data 
the  best  (within  1.7  percent  mean  error)  is 

C  /  Tr\  1/3 

-gf  =  2.72  +  0.263  (3.19) 

where  0.263  is  a  small  correction  term. 

The  maximum  heat  transfer  density  that  corresponds  to  the  optimal  spacing  of  Fig.  3.3 
is  reported  in  Fig.  3.6.  The  dimensionless  group  used  on  the  abscissa  of  Fig.  3.6  is  result 
of  the  theory:  it  is  easy  to  show  that  when  equation  (3.19)  is  substituted  without  the  0.263 
correction  into  equation  (3.4)  or  equation  (3.8),  the  group  qjRa]^'^  becomes  a  function  of  only 
the  group  {Hj Dyl^Ra^^^^  which  appeared  also  on  the  abscissa  of  Fig.  3.3.  The  numerical 
data,  however,  are  correlated  better  if  we  plot  q  instead  of  qfRa]^'^  on  the  ordinate.  The 
data  shown  in  Fig.  3.6  are  reproduced  within  1.7  percent  by  the  correlation 


qmax 


0.448 


1/3 


-1.6 


(3.20) 


3.4  Experimental  Results 

Figure  3.7  shows  the  main  features  of  the  arrays  used  in  the  experimental  phase  of  this 
study.  The  array  volume  was  fixed:  H  =  39.4mm,  L  =  142mm  and  W  =  44.5mm.  Three 
arrays  were  constructed  by  varying  from  5  to  3  the  number  of  horizontal  rows  in  the  H  xW 
cross-section.  The  number  of  cylinders  of  the  three  arrays  were  23,  14  and  8,  and  the 
respective  spacings  S/D  were  0.5,  1  and  2. 

The  longitudinal  conduction  and  loss  of  heat  through  the  cylinder  ends  were  minimized 
by  holding  each  cylinder  between  two  vertical  wooden  walls.  The  array  assembly  was  held 
inside  a  2m-high  enclosure  with  a  horizontal  cross-section  of  0.4m  x  0.4m.  The  bottom 
and  top  ends  of  the  enclosure  were  open  to  room  air.  To  reduce  the  effect  of  radiation,  the 
internal  and  external  surfaces  of  the  enclosure  were  covered  with  aluminum  foil.  By  using  the 
method  of  calculation  outlined  in  Morgan  (1975)  and  Sadeghipour  and  Asheghi  (1994),  we 
estimated  that  in  our  experiments  the  radiation  contribution  to  the  overall  heat  transfer  rate 
was  less  than  1.6  percent.  The  maximum  temperature  recorded  on  the  array  {T-wi,  Fig.  3.7) 
was  47.2°C. 

Each  cylinder  was  a  low  density  electrical  heater  consisting  of  a  helically  wound  heating 
element  (resistance  96  Q)  held  in  a  ceramic  insulator  filled  with  conductive  magnesium 
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Figure  3.7:  Details  of  the  test  section  of  the  experimental  apparatus. 


oxide.  The  outer  cover  was  polished  304  stainless  steel,  which  provided  adequate  rigidity 
and  resistance  to  oxidation.  The  heaters  were  connected  in  parallel,  and  powered  by  a 
variable  autotransformer  that  produced  a  voltage  between  0  and  140  V. 

The  temperatures  Tlu2  and  Too  were  measured  in  the  vertical  midplane  of  the  array, 
i.e.  half-way  between  the  two  wooden  walls.  Copper-constantan  type  T  thermocouples  were 
imbedded  in  1mm  hemispherical  depressions  machined  into  the  stainless  steel  sheath  of  the 
heater.  The  thermocouple  readings  were  referenced  to  a  mixture  of  crushed  ice  and  water. 
As  the  heat  transfer  rate  per  cylinder  was  distributed  uniformly  throughout  the  array,  the 
maximum  temperature  was  registered  on  the  trailing  row  of  cylinders  {T^x)-  The  temperature 
was  practically  uniform  (within  0.04°C)  around  the  cylinder  circumference:  we  measured  this 
variation  by  running  experiments  with  a  single  cylinder,  and  rotating  the  cylinder  to  change 
the  position  of  the  thermocouple  on  the  circumference. 

The  overall  thermal  conductance  q  was  calculated  after  measuring  the  power  dissipated 
in  all  the  heaters  (5),  the  maximum  temperature  (T^x)  and  the  room  temperature  at  the 
bottom  (inlet)  of  the  enclosure  (Too)-  The  q  values  were  calculated  using  T^i  in  place  of 
Tyj  in  equation  (3.18).  The  properties  of  air  were  evaluated  at  the  average  temperature 
{Tw  +  Tco)/2,  where  Tu,  =  {Tyjx  +  T^2)l2- 

The  uncertainties  associated  with  the  q  and  Ra£)  values  determined  in  this  manner  were 
2.5  percent  and,  respectively,  4.9  percent.  These  were  estimated  based  on  the  method  of 
Kline  and  McClintock  (1953),  and  the  following  inputs:  0.84  percent  uncertainty  in  T^ul,^^u2 
and  Too,  resulting  from  the  calibration  of  the  thermocouples;  0.5  and  1  percent  uncertainties 
in  the  measurement  of  voltage  and,  respectively,  current;  and  0.28  percent  for  Cp,  2  percent 
for  p,  and  1  percent  for  k. 

We  started  each  run  by  setting  the  voltage  and  current  for  the  resistance  heaters.  We 
then  waited  3  to  4  hours  while  monitoring  the  changes  in  voltage,  current,  T^^^^T^^  Too- 
We  ended  the  run  by  taking  final  readings  when  the  relative  change  in  voltage,  current  and 
temperature  readings  was  less  than  0.2  percent,  0.2  percent  and  0.06  percent  (determined 
by  repeating  10  hour-long  runs  at  the  same  spacing  and  Ran),  respectively.  These  relative 
changes  are  small  when  compared  with  the  uncertainties  in  the  respective  measurements. 

Figure  3.8  shows  a  summary  of  the  heat  transfer  experiments.  The  figure  shows  also 
the  curve  fitting  procedure  used  for  determining  S opt  ID,  i.e.  the  same  procedure  as  in  the 
numerical  part  of  the  study  (Fig.  3.5).  The  optimal  spacings  calculated  in  this  manner  are 
reported  in  Table  3.1.  Again,  the  best  practical  array  is  the  one  with  the  largest  overall 
thermal  conductance,  namely  the  case  n  =  14  (Figs.  3.7  and  3.8).  The  curve  fitting  of  the 
q  data  in  Fig.  3.8  is  used  only  as  a  test  of  the  accurancy  of  the  method  employed  in  the 
numerical  part  of  the  study. 


72 


We  used  these  experimental  results  to  verify  the  accuracy  of  the  numerical  method  de¬ 
scribed  in  the  preceding  section.  The  numerical  results  have  been  added  to  Table  3.1.  We 
performed  the  numerical  work  for  HjD  =  6.2  and  Pr  =  0.72,  which  correspond  to  the 
experimental  conditions.  First,  we  modelled  the  cylinders  as  surfaces  with  uniform  flux, 
and  found  that  the  maximum  temperature  occurs  at  the  top  of  the  topmost  cylinder.  The 
optimal  spacing  determined  numerically  agrees  within  17  percent  with  the  corresponding 
experimental  result. 

Table  3.1.  Comparison  between  the  optimal  spacings  determined  experimentally,  and 


the  corresponding  results  obtained  based  on  numerical  simulations  {HjD  = 
6.2,  Pr  =  0.72).  _ 


Ran 

Sopt/D 

200 

300 

350 

400 

Experiments 

1.44 

1.47 

1.45 

Numerical, 

uniform  flux  cylinders 

1.78 

1.68 

1.62 

1.61 

Numerical, 
isothermal  cylinders 

1.54 

1.48 

1.46 

1.44 

We  used  this  opportunity  to  see  if  the  type  of  cylinder  thermal  boundary  condition  has  an 
effect  on  the  optimal  spacing.  We  solved  again  the  experimental  configuration  numerically, 
this  time  modelling  the  cylinders  as  isothermal,  and  found  the  results  listed  in  the  bottom 
line  of  Table  3.1.  In  place  of  Ty,  in  equation  (3.18)  we  used  the  average  temperature  of  the 
topmost  cylinder.  It  is  clear  that  the  optimal  spacing  is  relatively  insensitive  to  whether  the 
array  is  uniform-fiux  or  isothermal.  This  conclusion  is  important,  because  it  means  that  the 
corelation  (3.19)  is  general. 

Additional  evidence  that  the  present  results  are  correct  is  provided  by  the  experimental 
study  of  Sparrow  and  Vemuri  (1985).  As  pointed  out  in  the  Introduction,  Sparrow  and 
Vemuri  maximized  the  total  heat  transfer  rate  from  an  isothermal  vertical  plate  [H  =  W  — 
7.62  cm)  with  horizontal  pin  fins  (Z)  =  0.635  cm).  Because  the  heat  transfer  maximum  was 
rather  flat.  Sparrow  and  Vemuri  concluded  that  the  optimal  spacing  occurs  when  the  total 
number  of  pins  fins  is  between  30  and  40.  The  flat  maximum  was  observed  at  Rayleigh 
numbers  (based  on  H)  in  the  range  4x10^  —  16x1 0^. 

Translated  into  the  terminology  employed  in  the  present  study.  Sparrow  and  Vemuri’s 
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conclusion  means  that  the  optimal  spacing  is  in  the  range  1.04  <  Soptj D  <  1.35,  when  the 
Rayleigh  number  is  in  the  range  230  <  Rao  <  920.  Their  SoptID  estimate  agrees  very 
well  (within  10  percent)  with  the  present  experimental  results  (Table  3.1),  especially  if  we 
recognize  that  (i)  the  H/ D  ratios  were  different  (12  for  Sparrow  and  Vemuri,  versus  6.2  in 
the  present  study),  (ii)  in  their  experiments  radiation  was  important,  and  (iii)  their  cylinders 
(pin  fins)  conducted  heat  longitudinally. 

3.5  Conclusion 

In  this  chapter  we  developed  fundamental  results  for  the  selection  of  the  spacing  between 
horizontal  cylinders  in  an  array  of  fixed  volume.  The  heat  transfer  is  by  laminar  natural 
convection.  The  optimal  spacing  reported  in  equation  (3.19)  corresponds  to  the  maximum 
thermal  conductance  between  the  entire  array  and  the  surrounding  fluid. 

More  important  fundamentally  is  the  conclusion  that  the  optimal  spacing  and  maximum 
thermal  conductance  can  be  expressed  in  compact  dimensionless  relations.  The  relevant  di¬ 
mensionless  groups  have  been  identified,  and  each  relation  contains  only  two  groups:  SoptID 
and  {Hf  for  optimal  spacing,  equation  (3.19),  and  ax  and  {HlDyl^Ra-j^''^  for 

maximum  thermal  conductance,  equation  (3.20). 

These  developments  were  possible  only  because  we  began  the  study  with  a  purely  theoret¬ 
ical  look  at  the  problem  of  predicting  the  optimal  spacing.  It  was  only  after  the  theory  that 
we  resorted  to  empiricism  (numerical  and  experimental),  which  was  necessary  in  order  to 
improve  the  accuracy  of  the  theoretical  expressions  found  for  SoptID  and  qmax-  Furthermore, 
the  dimensionless  groups  and  analytical  expressions  revealed  by  the  theory  had  the  effect 
of  minimizing  the  amount  of  numerical  and  experimental  information  needed  for  developing 
the  recommended  results,  equations  (3.19)  and  (3.20). 

3.6  Notation 

D  diameter 

g  gravitational  acceleration 

H  height  of  array.  Fig.  3.1 

k  thermal  conductivity 

L  cylinder  length.  Fig.  3.1 

n  number  of  cylinders  in  the  volume  HLW 

Tic  number  of  half-cylinder  in  the  computational  domain 
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Nujo  single  cylinder  overall  Nusselt  number 

p  pressure 

P  dimensionless  pressure 

Pr  Prandtl  number,  v j a 

qi  heat  transfer  from  a  single  cylinder 

q  heat  transfer  from  the  entire  array 

q  dimensionless  heat  transfer  volumetric  density 

Rud  Rayleigh  number,  g^D^{Tu,  —  T^)l{av) 

Ran  Rayleigh  number,  g/3H^{Tyj  —  Too)l{oii') 

S  cylinder-to-cylinder  spacing 

S  volume  averaged  spacing 

T  temperature 

Tw  cylinder  surface  temperature 

Too  coolant  inlet  temperature 

u,  V  velocity  components 

C/,  V  dimensionless  velocity  components 

U  mean  velocity  through  ^-wide  channel 

W  frontal  width  of  array.  Fig.  3.1 

X,?/  horizontal  and  vertical  coordinates 

X,  Y  dimensionless  coordinates 

a.  thermal  diffusivity 

coefficient  of  volumetric  thermal  expansion 
9  dimensionless  temperature 

p  density 
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Chapter  4 


THE  OPTIMAL  SPACING  FOR  CYLINDERS  IN  CROSS-FLOW 

FORCED  CONVECTION 


4.1  Introduction 

In  this  chapter  we  focus  on  new  fundamental  aspect  of  the  heat  transfer  performance  of 
a  bundle  of  parallel  cylinders  with  cross-flow  forced  convection,  namely,  the  maximization 
of  the  thermal  contact  between  the  bundle  and  the  fluid,  when  the  volume  occupied  by  the 
bundle  is  fixed.  In  the  experiments  described  by  Jubran  et  al.  (1993)  we  have  seen  empirical 
evidence  that  the  total  heat  transfer  rate  is  maximum  when  the  cylinder- to- cylinder  spacing 
S  has  a  certain  value.  This  finding  is  important  because  it  has  been  overlooked  for  decades, 
while  forced  convection  from  cylinders  in  cross-flow  grew  into  one  of  the  most  researched 
topics  in  heat  transfer. 

Jubran  et  al.  (1993)  did  not  offer  any  theoretical  explanation  or  way  of  predicting  the 
optimal  spacing.  The  objectives  in  this  chapter  are  two: 

(a)  To  show  that  the  optimal  SfD  ratio  can  be  predicted  based  on  a  simple  theory  (Bejan 
and  Sciubba,  1992),  and  that  contrary  to  Jubran  et  al.  ’s  conclusion,  the  SoptiD  ratio  is  not 
a  constant. 

(b)  To  also  show  that,  if  the  theory  (a)  is  known,  all  the  empirical  information  necessary 
for  predicting  SoptjD  accurately  is  already  available  in  the  large  volume  of  heat  and  fluid 
flow  data  published  for  cylinders  in  cross-flow. 

4.2  The  Method  of  Intersecting  the  Asymptotes 

It  is  useful  to  begin  with  the  simplest  analysis  that  (i)  proves  the  existence  of  an  optimal 
spacing,  and  (ii)  reveals  the  proper  trends  and  dimensionless  groups.  Consider  the  bundle 
shown  in  Fig.  4.1.  The  cylinders  occupy  the  fixed  volume  H  x  L  x  W,  where  H  is  aligned 
with  the  flow  direction.  The  cylinder  diameter  D,  the  pressure  drop  across  the  bundle  AP 
and  the  upstream  temperature  Too  are  also  fixed.  The  cylinder  temperature  is  of  order  T^,; 
this  order  of  magnitude  characterizes  all  the  cylinders  in  the  bundle.  We  are  interested 
in  maximizing  the  total  heat  transfer  q  between  the  bundle  and  the  surrounding  fluid,  by 
selecting  the  cylinder-to-cylinder  spacing  S,  or  the  number  of  cylinders  in  the  bundle. 

Consider  first  the  limit  where  the  spacing  S  is  sufficiently  large  that  each  cylinder  acts 
as  if  it  is  alone  in  its  own  cross-flow  of  free-stream  velocity  Uoo-  The  total  heat  transfer  rate 
experienced  by  the  bundle  is  g  =  ngi,  where  qi  is  the  heat  transfer  associated  with  a  single 
cylinder. 
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Figure  4.1:  Fixed  volume  H  xLxW  containing  a  bundle  of  parallel  cylinders  perpendicular 
to  a  free  stream. 
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(4.1) 


9i  =  ^NuirDL{Tw  -  Too) 
and  n  is  the  total  number  of  cylinders, 


(5  +  D)2cos30°  ^  ^ 

We  are  assuming  that  W  is  considerably  greater  than  (S  +  D).  In  the  range  0.7  <  Pr  <  500 
and  40  <  UooDjv  <  1000,  the  average  Nusselt  number  is  given  by  the  correlation  (Zukauskas, 
1987b) 


iVu  =  0.52Pr°-^^i?e^/^  (4.3) 

The  free  stream  velocity  Uoo  is  not  given.  It  is  determined  by  the  force  balance  on  the 
entire  bundle,  AP  •  WL  =  nPi,  where  Pi  is  the  drag  force  experienced  by  one  cylinder, 
Pj  =  CdDL{pU^I2).  The  drag  coefficient  varies  from  2  to  1  in  the  Re  range  40  -  1000, 
therefore  in  this  order  of  magnitude  analysis  it  is  sufficient  to  use  Cd  ~  1-5.  The  force 
balance  yields  Uoo,  or  Re, 


/  n  \  n  TT 

1.1(5 +  l)(ij  ,  with  S=-,  H  =  ^  (4.4) 

where  P  is  the  dimensionless  pressure  group  identified  by  Bhattacharjee  and  Grosshandler 
(1988) 


P  _  AP  • 

fJLU 

Combining  equations  (4.1)  -  (4.6)  we  find  that  the  total  heat  transfer  rate  behaves  as 


(4.5) 


W 

qiargeS  —  2A;T  — (T^ 


Too) 


^3/4pl/4p^0.37 

(P  +  1)3/2 


(4.6) 


Consider  now  the  opposite  extreme  when  the  cylinders  almost  touch,  and  the  flow  is 
almost  cut  off.  In  this  limit  the  temperature  of  the  coolant  that  exits  slowly  through  the 
right  end  of  the  bundle  (the  plane  L  x  W)  is  essentially  the  same  as  the  cylinder  temperature 
Pu,.  The  heat  transfer  from  the  bundle  to  the  coolant  is  equal  to  the  enthalpy  gained  by  the 
coolant,  q  =  mcp{T^  —  Too)  ,  where  m  is  the  mass  flowrate  through  the  L  xW  plane. 

To  obtain  an  order-of-magnitude  estimate  for  the  flowrate,  we  note  that  m  is  composed 
of  several  streams  [total  number  rit  =  W/{S  +  D)],  each  with  a  cross-sectional  area  S'  x  P  in 
the  plane  of  one  row  of  cylinder  axes.  The  thickness  of  the  channel  traveled  by  each  stream 
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varies  between  a  minimum  value  {S)  at  the  row  level,  and  a  maximum  value  at  a  certain 
level  between  two  rows.  The  volume- averaged  thickness  of  one  channel  of  this  kind  is 


S  =  S  +  B  -  0.907^^  (4.7) 

however,  we  may  adjust  this  estimate  by  using  1  in  place  of  the  factor  0.907  to  account  for 
the  fact  that  the  flow  must  cease  when  the  cylinders  touch  (S  =  0): 


S  =  S 


S  +  2D 
S  +  D 


(4.8) 


The  mean  velocity  U  through  a  channel  of  length  H  and  cross-sectional  area  5T  can  be 
estimated  using  the  solution  for  Hagen-Poiseuille  flow  through  a  parallel  plate  channel  of 
spacing  S  and  length  if. 


S^AP 

l2fiH 


(4.9) 


The  mass  flowrate  through  one  channel  is  mi  =  pUSL.  There  are  Ut  channels,  therefore  the 
mass  flowrate  through  the  entire  bundle  is  m  =  miWI{S  +  D),  and  q  becomes 


1  W 

qsmallS  =  ~ 


PPrS^iS  +  2)3 
H{s  +  iy 


(4.10) 


The  two  asymptotic  trends  are  sketched  in  Fig.  4.2.  The  actual  (unknown)  curve  g(5), 
which  is  indicated  by  the  solid  line  in  the  figure,  has  a  maximum  where  the  spacing  S  is 
approximately  the  same  as  the  S  value  obtained  by  intersecting  the  two  asymptotes  (Bejan 
and  Sciubba,  1992).  The  Sopt  value  obtained  by  eliminating  q  between  equations  (4.6)  and 

(4.10)  is  given  implicitly  by 


,  (2 + sj'.ri-’  ^ 

"’"(l +5„-i)»/’  '  p3/upr^M 


(4.11) 


The  optimal  spacing  increases  with  the  length  of  the  bundle,  and  decreases  with  the 
applied  pressure  difference  and  the  Prandtl  number.  It  is  also  interesting  that  in  equation 

(4.11)  the  exponents  of  P  and  Pr  are  almost  the  same  (note  that  3/14  =  0.21).  This  means 
that  instead  of  the  product  p3/i4py,o.i8  approximately  II^/i^,  Fig.  4.3, 


,  (2  +  ^ 

°"(i  +  5-i)'/'  ■  nvH 


(4.12) 


where  the  11  group  is  defined  as  11  =  AP  ■  D'^Kfia)  =  PPr,  (Bejan,  1993).  The  intersection 
of  the  large-S"  and  small-5  regimes  provides  also  an  estimate  for  the  scale  of  the  maximum 
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Figure  4.2:  The  optimal  cylinder-to-cylinder  spacing  as  the  intersection  of  the  large-F  and 
small- 5"  asymptotes. 


thermal  conductance  between  bundle  and  coolant,  ql{Tw  —  Too)-  Figure  4.2  shows  that  the 
result  of  this  operation  will  always  overestimate  the  peak  value  of  the  actual  ql{Tw  —  Too) 
curve: 


kLW 

~w 


{T^ 


Too) 


>  1-85 


(4.13) 


4.3  Estimates  Based  on  Experimental  Correlations  for  Cross-Flow 

Now  that  we  have  theoretical  reasons  to  expect  an  optimal  spacing  and  a  scaling  of 
type  (4.13),  we  turn  our  attention  to  a  more  precise  analysis  that  is  based  on  published 
experimental  data.  A  very  large  effort  has  been  devoted  to  documenting  the  pressure  drop 
and  heat  transfer  coefficient  for  bundles  of  parallel  cylinders  in  cross-flow  (e.g.  Zukauskas, 
1987b).  These  results  are  expressed  in  terms  of  the  Reynolds  number  based  on  the  maximum 
average  velocity  between  adjacent  cylinders,  Rcmax  =  VmaxDjv.  To  determine  the  relation 
between  the  applied  pressure  difference  and  Rcmax  we  use 


AP  =  nifx]^pylax  (4-14) 

with  X  =  1  for  the  equilateral  triangle  array,  and  ni  for  the  number  of  rows  in  the  flow 
direction,  nc  =  I  {H  -  1)/0.866(P  +  1).  Equation  (4.14)  becomes 


P  =  \rl>fR^L.  (^-15) 

The  friction  factor  /  is  available  in  graphic  form  (see  Fig.  6.20  in  Zukauskas,  1987b)  as  a  func¬ 
tion  of  Rcmax  and  cylinder-to-cylinder  spacing.  In  the  range  10  <  Rtmax  <  10^  and  0.25  < 
S  <  1.5,  the  /  curves  are  approximated  within  10  percent  by  /  =  {G\l Rtmax)  +  C'2,  where 
Cl  =  68,S“°'^®  and  C2  =  0.375“°'^^.  Equation  (4.15)  provides  the  function  Remax{P,S)  that 
accounts  for  the  flow  part  of  the  problem. 

When  the  total  heat  transfer  q  is  divided  equally  among  the  cylinders  the  peak  tempera¬ 
ture  {Tmax)  belongs  to  the  cylinders  positioned  in  the  last  row.  The  maximum  temperature 
difference  {Tmax  —  Too)  can  be  calculated  in  two  steps  (Knight  et  al,  1991), 


Tmax  Tout 


{qln)D 


and  Tout  Too  — 


wDLkNu  mcp 

where  Tout  is  the  bulk  temperature  of  the  fluid  leaving  the  bundle,  m  =  pVm 
total  mass  flowrate,  and  ut  =^Wl{S  D).  Eliminating  Tout  we  obtain 


(4.16) 
rSLut  is  the 


T  = 


kLW 

qD 


{Trr. 


-Too) 


5+1  ^  5+1 

-KTliNu  SPrRemax 


(4.17) 
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where  T  is  the  dimensionless  maximum  temperature  difference.  Finally,  we  assume  that 
there  are  many  cylinders  in  the  longitudinal  flow  direction  (n^  >  16),  such  that  the  Nusselt 
number  for  the  last  row  is  the  same  as  the  value  averaged  over  the  entire  bundle.  The 
numerical  results  reported  below  were  developed  using  Table  6.3  of  Zukauskas  (1987b). 

In  conclusion,  the  result  of  combining  equations  (4.15)  and  (4.17)  is  a  maximum  tempera¬ 
ture  difference  function  f  that  depends  on  P,  H,  Pr  and  S.  The  f  function  can  be  minimized 
numerically  with  respect  to  the  spacing  S,  to  obtain  the  optimal  spacing  Sopt{P,  H,  Pr).  Fig¬ 
ure  4.3  shows  the  Sopt  results  obtained  for  a  total  of  36  P,  H  and  Pr  combinations  in  the 
range  lO"*  <  P  <  10*,  25  <  H  <  200,  and  0.72  <  Pr  <  50.  In  conclusion,  the  estimate 
obtained  in  the  preceding  section  is  adequate  as  intended  (i.e.  in  an  order  of  magnitude 
sense),  and  the  predicted  trends  are  also  correct. 

The  scatter  of  the  points  plotted  in  Fig.  4.3  suggests  that  the  use  of  the  combined  group 
jji/2/jp/i4  abscissa  may  not  be  the  best  way  to  correlate  the  empirical  data.  After 

examining  individually  the  effects  of  H,P  and  Pr  on  Sopt,  we  found  that  the  following 
correlation  approximates  the  data  with  a  standard  deviation  of  5.6  percent: 


ffO.52 

Sopt  =  1.59^ - 

^  pO.13  ppO.24 


(4,18) 


The  minimum  last-row  (i.e.  peak)  temperature  occurs  when  the  spacing  is  optimal,  Tmin  = 
T{Sopt)-  The  Tmin  values  for  the  36  cases  documented  in  Fig.  4.3  are  correlated  with  a 
standard  deviation  of  16  percent  by  the  expression 


3.33 


Tmin  -  ^0.45pr0.64 

In  both  equations  (4.13)  and  (4.19),  the  denominator  on  the  right-hand  side  is  approximated 
well  by  This  allows  us  to  rewrite  equation  (4.19)  as  the  maximum  power  density 

installed  in  the  fixed  volume  HLW, 


q 

HLW 


0  311^^^^^^’^°^ 


HD 


(4.20) 


We  note  that  11  is  proportional  to  therefore  the  maximum  power  density  does  not  depend 
on  the  cylinder  diameter.  Important  is  to  set  the  spacing  S  in  proportion  to  D,  c.f.  equation 
(4.18). 


4.4  Concluding  Remarks 

It  is  instructive  to  compare  the  optimal  spacing  reported  by  Jubran  et  al.  (1993)  with  the 
data  plotted  in  Fig.  4.3.  In  the  present  notation,  Jubran  et  a/.’s  conclusion  reads  Sopt  =  1-89, 
and  was  reached  during  experiments  with  air  flow  through  arrays  of  parallel  pin  fins  ( D  = 
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6.35  mm,  L  =  60  mm)  attached  to  a  base  plate  {H  =  300  mm,  W  =  175mm).  A  direct 
comparison  with  the  present  results  is  not  possible  because  in  Jubran  et  a/.’s  experiments 
S  was  varied  at  constant  air  mass  flowrate  (or  constant  upstream  velocity,  Uin),  instead  of 
constant  AP.  The  experiments  covered  the  Uin  range  2.2  -  lOm/s. 

An  order  of  magnitude  comparisons  can  be  made,  however,  by  noting  that  when  5  ~  2 
the  array  is  dense  and  long  enough  such  that  the  overall  AP  scale  is  AP  ~  The  array 

appears  “dense”  from  the  point  of  view  of  the  approaching  stream,  to  which  the  array  looks 
opaque  (without  clear  openings).  The  stream  has  two  choices,  to  penetrate  the  array,  or  to 
flow  around  it  (Morega  et  ah,  1994).  It  does  both,  with  less  penetration  as  S  decreases. 

If  we  use  this  approximation  we  find  that  Jubran  et  a/.’s  experiments  cover  a  range 
characterized  by  4.2  x  10®  <  P  <  9  x  10®,  H  —  47.3  and  Pr  —  0.72.  Their  conclusion  that 
Sopt  =  1.89  is  shown  as  a  horizontal  bar  in  Fig.  4.3,  and  agrees  very  well  with  the  Sopt  values 
determined  by  minimizing  the  f  expression  of  equation  (4.17).  This  agreement  is  remarkable 
especially  in  view  of  the  fact  that  in  Jubran  et  a/.’s  experiments  the  cylinders  were  pin  fins 
attached  to  a  base  surface,  i.e.  the  cylinders  were  not  long  relative  to  their  diameter,  their 
temperature  varied  along  their  length,  and  the  time-averaged  flow  was  not  two-dimensional. 

The  same  AP  approximation  can  be  used  to  convert  the  empirical  correlations  (4.18) 
and  (4.19)  into  order-of-magnitude  formulas  for  bundles  exposed  to  a  free  stream  of  velocity 


4.5 

PePpPrO-64 


(4.21) 


where  Rdn  =  UinD/u,  which  is  the  same  as  Rein  ~  (2P)^/^.  The  correlations  (4.21),  cover 
the  range  140  <  Pe,n  <  14000,  25  <  H  <  200,  and  0.72  <  Pr  <  50. 

The  relative  progress  made  in  this  chapter  is  illustrated  in  Fig.  4.3.  Known  until  now 
was  the  single  point  reported  by  Jubran  et  al.  (1993):  the  point  (horizontal  bar)  shown 
in  Fig.  4.3,  however,  is  more  general  because  the  cartesian  coordinates  are  dimensionless. 
These  coordinates  are  the  product  of  this  chapter.  The  other  product  is  the  complete  curve, 
equation  (4.18),  which  is  based  on  experimental  heat  and  fluid  flow  data,  and  covers  all  the 
possible  flow  conditions. 


4.5  Notation 


Ci,2  functions  of  S 
Cd  drag  coefficient 
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D 

f 

F, 

h 

H 

H 

k 

L 

m 

n 

ni 


nt 


Nu 

P 


Q 


9i 


cylinder  diameter 

friction  factor 

drag  force  per  cylinder 

average  heat  transfer  coefficient 

bundle  flow  length 

dimensionless  length,  H/D 

coolant  thermal  conductivity 

cylinder  length 

total  mass  flowrate 

total  number  of  cylinders 

number  of  rows  in  the  longitudinal  direction 

number  of  rows  in  the  transversal  direction 

average  Nusselt  number,  hD  jk 

pressure  drop  number,  AP  •  ((nv) 

total  heat  transfer  rate 

heat  transfer  rate  per  cylinder 

Reynolds  numbers,  UooD/v,  UinDjy,  VmaxD/v 

cylinder-to-cylinder  spacing 

dimensionless  spacing,  S/D 

average  spacing 

dimensionless  maximum  temperature  difference 
coolant  outlet  temperature 
cylinder  temperature  scale 
coolant  inlet  temperature 
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U  mean  velocity 

Uin  inlet  free-stream  velocity 

Uoo  free-stream  velocity  around  one  cylinder 

Vmax  maximum  average  velocity 

W  bundle  width 

a  coolant  thermal  diffusivity 

AP  pressure  difference 

fjL  viscosity 

V  kinematic  viscosity 

n  pressure  drop  number,  AP  ■  Z)^/()Ua:) 

p  coolant  density 

X  correction  factor 
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Chapter  5 


THE  OPTIMAL  SIZES  OF  BODIES  WITH  EXTERNAL  FORCED 
CONVECTION  HEAT  TRANSFER 


5.1  Introduction 

External  forced  convection  is  one  of  the  most  common  heat  transfer  modes  that  occur  in 
the  heat  exchangers  used  in  the  power  and  refrigeration  industry  (e.g.  plate  and  cylindrical  fins, 
tubes  in  cross-flow).  It  is  also  well  known  that  each  heat  transfer  process  is  accompanied  by 
entropy  generation,  which  a  quantity  proportional  to  the  corresponding  useful  work  (exergy) 
destroyed  by  the  power  or  refrigeration  plant.  The  ways  to  minimize  the  generation  of  entropy  in 
such  processes  have  become  a  fundamental  topic  in  contemporary  heat  transfer. 

The  objective  of  this  chapter  is  to  show  how  to  optimize  the  size  of  a  body  engaged  in 
external  forced  convection,  when  its  heat  transfer  duty  is  specified.  Two  optimization 
approaches  are  used.  The  first  is  based  on  minimizing  the  total  rate  of  entropy  generation  due  to 
fluid  friction  and  imperfect  thermal  contact.  In  the  second  approach,  the  minimized  quantity  is 
the  total  cost  associated  with  the  material  of  the  body  and  the  imperfect  heat  transfer. 


5.2  Optimal  Sizes  for  Minimum  Entropy  Generation 

The  rate  of  entropy  generation  associated  with  external  forced  convection  is,  in  general 
(Bejan,  1982,  1988), 


q(Tw-Tj 

_  O  * 


FpU 

T^ 


(5.1) 


In  writing  equation  (5.1)  it  is  assumed  that  the  surface- averaged  temperature  difference  between 
the  body  and  the  free  stream  (T^  —  T„)  is  at  least  one  order  of  magnitude  smaller  than  the 


absolute  temperature  Too.  The  two  terms  on  the  right-hand  side  of  equation  (5.1)  represent,  in 
order,  the  contribution  due  to  heat  transfer  across  a  finite  temperature  difference,  and  the 
contribution  due  to  fluid  friction.  The  first  applications  of  equation  (5.1)  were  to  the  flat  plate 
with  laminar  and  turbulent  boundary  layer  flow  (Bejan,  1982,  Example  5.3  and  Problem  5.4)  and 
the  sphere  (Bejan,  1982,  Problem  5.3),  where  it  was  shown  that  the  body  dimension  can  be 
selected  such  that  is  minimized.  The  cylinder  in  cross-flow  and  the  flat  plate  with  turbulent 

boundary  layer  flow  were  optimized  by  Poulikakos  and  Johnson  (1989).  We  begin  the 
correlation  work  by  extending  the  existing  results  to  wider  Reynolds  number  ranges. 

5.2.1  The  smooth  cylinder  in  cross-flow.  Consider  the  heat  transfer  and  friction 
between  a  long  cylinder  of  diameter  D  and  length  W,  and  a  uniform  cross-flow  (Uoo,  T«,).  The 
total  heat  transfer  rate  (the  heat  transfer  duty)  is  fixed,  q  =  h7cDW(T^  -T^),  where  T^  is  the 

perimeter  averaged  temperature.  The  heat  transfer  coefficient  h  can  be  estimated  based  on 
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Churchill  and  Bernstein's  (1977)  correlation,  which  is  valid  for  all  values  of  Reo  and  Pr, 
provided  ReoPr  >  0.2.  The  total  drag  force  is  Fd  =  CdDW  pUif  /  2,  for  which  the  relationship 
between  the  drag  coefficient  Cd  and  Reo  is  available  in  graphical  form  (e.g.  Bejan,  1993,  p. 
266).  We  digitized  and  curve-fitted  the  CdCRcd)  data  in  the  Reo  range  0.1  -  10  5,  and 
substituted  this  information  into  the  second  term  on  the  right  side  of  equation  (5.1), 


/,f--^4RepC„(Re„)  (5.2) 

HUi/T.  itNu(ReB,Pr)  2  “  “ 

where  Nu  =  hD  /  k  and  B  is  the  nondimensional  heat  transfer  duty  parameter  of  the  cylinder. 


B=  ig  (5.3) 

We  minimized  the  expression  (5.2)  numerically  by  choosing  the  optimal  diameter 
(dimensionless  Reo,  opt)  for  a  given  (B,  Pr)  pair.  The  B  definition  (5.3)  was  chosen  in  such  a 
way  that  the  results  can  be  projected  as  a  single  curve  Reo,  opt  (B)  in  Fig.  5.1,  if  the  Prandtl 
number  is  0.7  or  greater. 

5.2.2  The  smooth  sphere.  The  same  optimization  method  can  be  applied  to  a  sphere  (D) 
with  heat  transfer  to  an  external  flow  (Uoo,  Too).  For  conciseness,  we  show  only  the  beginning 
and  the  end  of  the  analysis.  The  overall  Nusselt  number  can  be  estimated  based  on  Whitaker's 
(1972)  correlation,  which  has  been  tested  for  0.71  <  Pr  <  380  and  3.5  <  Reo  <  7.6  x  10^.  The 
drag  force  is  Fp  =  (tt  D’^  /  4)  p  /  2  where  CD(ReD)  is  the  drag  coefficient  (Bejan,  1993). 
After  making  this  substitution  for  Fd  in  equation  (5.1)  we  were  able  to  write 


'gen 


B^Pr^ 


+ 


v2  p  /  T^  %  Reo  Nu  (Rep,  Pr) 
where,  Nu  =  hD  /  k,  and  the  given  heat  transfer  rate  is 


^  Rep  Cd  (Rej^)  (5.4) 

q  =  h  tcD^  (T^  -T^)  .  We  minimized  this 


expression  numerically  with  respect  to  Rcd,  and  plotted  the  resulting  ReD,  opt  values  versus  Bg  in 
Fig.  5.1.  Once  again,  Bg  was  chosen  such  that  the  results  are  represented  by  a  single  curve  when 
Pr  >  0.7.  The  sphere  heat  transfer  duty  parameter  Bg  is  defined  by 


B,=— - (5.5) 

v(knT„Pr''’)‘'^ 

5.2.3  The  flat  plate  in  parallel  flow.  Projected  on  Fig.  5.1  are  also  the  results  known 
for  the  flat  plate  of  length  L,  with  the  heat  flux  q"  on  both  sides.  The  total  heat  transfer  rate  is 
fixed,  q  =  2LWq",  where  W  is  the  plate  width.  Using  the  present  nondimensionalization 
approach,  the  results  known  for  laminar  boundary  layer  flow  (ReL  <  5  x  10^)  can  be  summarized 
here  as 
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@2 


=  0.554  B2, 


(5.6) 


Re 


L,opt 


S  .  /W 

*^gen,mm^  9  -d 


For  the  optimization  of  the  flat  plate  in  the  turbulent  boundary  layer  regime,  a  wider  ReL 
range  (5  x  10^  -  10^)  is  covered  by  the  correlation  (e.g.  Bejan,  1993,  p.  261): 
q"L/k(T^-T„)  =  0.037Pr^/^(ReJ'^-23,550).  To  estimate  the  drag  force  Fd  =  2LW  x  we 
used  X  =  0.037  p  Rcl^^^  .  The  optimal  Rol  that  minimizes  the  resulting  Sgg^  expression  is 


(Fig.  5.1) 


Rei., ^  =  (13.516+23,550)*'“ 


(5.7) 


5.2.4  Effective  length,  and  universal  duty  parameter.  The  main  contribution  of  the 
dimensionless  chart  constructed  in  Fig.  5.1  is  to  show  that  the  optimal  size  of  the  object 
(Rcd,  opt.  ReL,  opt)  increases  monotonically  with  its  respective  duty  parameter  (B,  Bg).  The 
challenge  then  is  to  identify  a  more  general  set  of  coordinates  in  which  the  three  curves  of  Fig. 
5.1  fall  on  the  same  curve. 

A  step  in  this  direction  is  presented  in  Fig.  5.2.  The  parameter  that  is  used  now  on  the 
ordinate  for  the  flat  plate  is  the  Reynolds  number  based  on  the  boundary  layer  thickness  at  the 
trailing  edge,  Res^opt  =  U  5opt/v.  We  made  this  choice  because  5  is  the  transversal  length 
scale  of  the  flow,  i.e.  the  same  type  of  flow  length  scale  as  the  thickness  of  the  wake  (D)  behind 
the  cylinder  and  the  sphere.  For  the  laminar  regime  on  the  flat  plate  we  used  the  estimate 
5  =  4.92  L  Refj^^^,  so  that  equation  (5.6)  was  replaced  by  Reg,  opt  =  3.66B.  For  the  turbulent 
regime  we  used  5  =  0.37  L  Rbl^^^  and  replaced  equation  (5.7)  with  Res,  opt  =  0-37  (13.51B  + 
23,550). 

The  abscissa  of  Fig.  5.2  shows  the  duty  parameter  B  of  equation  (5.3),  which  was  defined 
in  the  same  way  for  the  cylinder  and  the  plate,  and  the  new  parameter  B*  for  the  sphere.  The 
latter  is  based  on  the  observation  that  the  B  and  Bg  definitions  (5.3)  and  (5.5)  differ  with  respect 
to  the  way  in  which  the  heat  transfer  duty  is  specified,  q  /  W  as  opposed  to  q.  We  can  construct  a 
sphere  duty  parameter  of  type  B,  by  regarding  the  spherical  surface  TtD^  as  equal  (in  an  order  of 
magnitude  sense)  to  the  lateral  surface  of  a  short  cylinder  (Do  =  W)  in  cross-flow  and  with 
insulated  ends.  From  the  area  conservation  statement  TtD^  =  ttDcW  we  learn  that  Dc  =  D,  in 


other  words,  the  diameter  of  the  equivalent  short  cylinder  is  the  same  as  that  of  the  actual  sphere. 
By  substituting  D  in  place  of  W  in  the  B  definition  (5.3),  we  arrive  at  the  B-t5rpe  duty  parameter 


for  the  sphere 

q/D  B, 

(k  pT„.  ReD.opt 


(5.8) 


Figure  5.2  shows  that  the  consistent  use  of  the  transversal  length  scale  and  the  duty 
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Figure  5.2  The  optimal  sizes  of  Fig.  5.1  correlated  in  terms  of  a  universal  duty  parameter  (B, 
B*)  and  the  Reynolds  number  based  on  the  transversal  length  scale  of  the  flow. 


parameter  of  type  B  brings  the  curves  much  closer  than  in  the  original  plot  of  Fig.  5.1.  This 
discovery  suggests  that  the  unknown  (and,  presumably,  more  complicated)  curve  that  represents 
a  body  of  another  shape  will  fall  not  far  from  the  curves  of  Fig.  5.2,  provided  the  ordinate  and 
abscissa  parameters  are  defined  in  a  manner  consistent  with  those  of  Fig.  5.2. 

5.3  Optimal  Sizes  for  Minimum  Cost 

5.3.1  When  heat  transfer  losses  dominate.  To  see  the  reasons  for  the  alternative 
approach  followed  in  this  section,  we  begin  with  a  numerical  example.  The  question  pursued  in 
this  example  is  how  the  optimal  sizes  of  Fig.  5.1  compare  with  the  usual  sizes  that  occur  in  real- 
life  applications.  It  is  sufficient  to  use  a  textbook  example  like  the  one  given  on  p.  415  in 
Incropera  and  DeWitt  (1990).  A  cylinder  with  D  =  1.27  cm,  W  =  9.4  cm  and  q  =  39.1  W  is 
cooled  in  a  cross-flow  of  air  with  Uoo  =10  m/s  and  Too  =  300  K.  If  we  evaluate  the  air  properties 
at  the  fUm  temperature  335  K,  the  duty  parameter  calculated  with  equation  (5.3)  is  B  =  3400.  We 
place  this  value  on  the  abscissa  of  Fig.  1,  and  read  ReD_opt  =  3.5  x  10“^,  from  which  we  deduce 
that  Dopj  =  6.6  cm. 

In  conclusion,  in  this  example  the  actual  diameter  (1.27  cm)  is  roughly  one  order  of 
magnitude  smaller  than  the  optimal  diameter  for  minimum  total  entropy  generation.  When  D  < 
Dopt  the  total  entropy  generation  rate  is  dominated  by  the  contribution  made  by  heat  transfer. 
The  question  then  is  how  to  determine  the  size  of  an  object  with  forced  convection  heat  transfer 
when  its  irreversibility  is  dominated  by  heat  transfer  irreversibility.  The  way  to  reduce  the  heat 
transfer  irreversibility  when  the  heat  transfer  duty  is  fixed  is  by  increasing  the  size  (or  wetted 
surface)  of  the  immersed  object.  This  increase  in  size  is  always  accompanied  by  an  increase  in 
the  cost  of  the  object  itself.  We  arrive  in  this  way  at  a  trade-off  between  saving  exergy  and 
spending  more  on  a  larger  object. 

In  this  section  we  examine  this  trade-off  by  putting  both  the  destroyed  exergy  and  the 
object  on  a  cost  basis.  According  to  equation  (5.1),  the  entropy  generation  rate  due  to  heat 
transfer  is  q  (T^  -  T„)  /  t£  We  multiply  this  by  the  absolute  temperature  Too  and  obtain  the 

corresponding  rate  of  exergy  destruction  in  the  installation  that  employs  the  object  as  a  heat 
transfer  device,  q  (T^  -  T„)  /  T^.  The  exergy  destroyed  during  the  lifetime  of  operation  to  is 
to  9  {^w  “  Too)  /  To,.  The  cost  associated  with  this  quantity  is  c^  t^,  q  (T^  -  T,,)  /  T,,,  where  Ce($/J) 
is  the  unit  cost  of  exergy  (e.g.  electric  energy  from  the  network). 

5.3.2  The  smooth  cylinder  in  cross-flow.  If  the  object  (cylinder,  in  this  case)  is  solid, 
its  cost  is  proportional  to  its  volume,  Cni(7t/4)  D^W,  where  Cm($  /  m^)  is  the  unit  cost  of  the 
cylinder  material.  This  unit  cost  may  refer  to  the  purchase  of  the  material  if  the  material  is 
expensive,  or  to  its  weight  if  the  entire  installation  is  to  be  part  of  an  aerospace  application.  In 
summary,  the  two  cost  items  can  be  linked  into  an  additive  cost  formula. 
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which  can  be  nondimensionlized  as 


C/W  _  Pr 


Cm(v/Uj^  ^Nu  4 
The  dimensionless  heat  transfer  duty  parameter  is  now 


+  -?Re2 


D 


(5.10) 


F  = 


q/W 


Ceto 


1/2 


(5.11) 


V  /  u„-  1  ^  l'  T  Pt-  1/3 

The  Nusselt  number  Nu  (Reo,  Pr)  is  given  by  the  Churchill  and  Bernstein  (1977)  correlation. 
The  optimal  Reynolds  number  for  which  the  total  cost  of  equation  (5.10)  is  minimized  emerges 
as  a  single  curve  of  Rep^  opt  versus  F  when  Pr  >  0.7.  This  curve  is  shown  in  Fig.  5.3. 

5.3.3  The  smooth  sphere.  The  minimization  of  the  cost  associated  with  a  solid  sphere 
can  be  performed  in  a  similar  way.  The  second  term  on  the  right  side  of  equation  (5.9)  is 
replaced  by  Cm  (7t/6)D3,  so  that  the  dimensionless  total  cost  formula  becomes 


where 


C 

Cm  ^  ^  ^  7t  Rej)  Nu 


+  fRe 

o 


3 

D 


(5.12) 


1  v  J 

iC„kT.Pr'»J 

1/2 


(5.13) 


By  using  the  Nu  (Reo,  Pr)  correlation  given  by  Whitaker  (1972),  it  is  possible  to  minimize  the 
expression  (5.12)  with  respect  to  Reo-  The  results  can  be  plotted  as  a  single  curve  if  Pr  >  0.7, 
Fig.  5.3. 

5.3.4  The  flat  plate.  If  the  boundary  layer  flow  is  laminar  and  the  heat  flux  q"  is 
uniform  on  both  sides  of  the  plate,  the  two-term  cost  formula  reduces  to 


=  0.736  f2  Ref  1/2 +  ReL  (Pr>0.5)  (5.14) 

where  t  is  the  plate  thickness,  t  «  L,  and 


r  q/W  f  Cet,  W/2 
(tv/uj^^  U^kT^Pri/3j 
The  optimal  plate  length  that  minimizes  the  total  cost  is 


(5.15) 


ReL,opt  =  0-514F4^  (5.16) 

For  cases  in  which  the  flow  is  turbulent,  we  use  the  correlation  listed  under  equation  (5.6) 
and  obtain  a  relation  similar  to  equation  (5.14),  except  that  the  first  term  on  the  right  side  is 
replaced  by  13.51  F^/(ReL^^- 23,550).  The  end  result  is  that  the  optimal  plate  length  for 
minimum  cost  is  described  implicitly  by 


F  =  0.304  (Rcl*"-  23,550  ) 


(5.17) 


where  the  duty  parameter  F  continues  to  be  defined  as  in  equation  (5.15).  Both  equations  (5.16) 
and  (5.17)  have  been  added  to  Fig.  5.3. 

5.3.5  Correlation  of  the  optimal  sizes  of  the  cylinder,  sphere  and  plate.  The  results 
assembled  in  Fig.  5.3  look  remarkably  similar  to  those  in  Fig.  5.1.  These  similarities  are  due  to 
the  fact  that  the  new  dimensionless  duty  parameters  (F,  Fg)  were  carefully  chosen  so  that  they 
have  the  same  structure  as  the  abscissa  parameters  of  Fig.  5.1.  This  means  that  we  should  be 
able  to  correlate  the  three  curves  of  Fig.  5.3  into  a  universal  curve,  by  following  the  method  of 
section  5.2.4,  which  led  to  Fig.  5.2. 

First,  in  place  of  the  sphere  duty  parameter  Fg  we  use  the  equivalent  parameter  of  type  F 
[equation  (5.11)]  by  imagining  a  short  cylinder  (D  -  W)  that  has  a  cylindrical  area  equal  to  that 
of  the  sphere, 


q/D  f  c,t,  11/2  _ 
'^*"v/U^[c^kT„Pri/3j  -  Reo 


(5.18) 


The  relation  between  the  abscissas  of  Figs.  5.3  and  5.4  for  the  sphere  is  F*  =  Fs  /  ReD,opt- 
Second,  for  the  plate  we  base  the  Reynolds  number  on  the  transversal  length  scale  of  the  flow, 
which  is  5  =  4.92  L  Ref^^  in  the  laminar  regime,  and  5  =  0.37  L  Ref^^^  in  the  turbulent  regime. 
The  ordinate  for  the  plate  changes  to  Re§,  opt  =  UooSopt/ v,  and  equations  (5.16)  and  (5.17)  are 
replaced  by  Reg^  =  3.53  F  and,  respectively,  F  =  0.0173  Re^^^pj  -  5206  Regy^p^.  Figure  5.4 

shows  once  again  that  this  method  of  correlation  (transversal  length  scale  on  the  ordinate  and 
cylinder-type  duty  parameter  on  the  abscissa)  works  very  well. 


5.4  Conclusions 

In  this  chapter  we  showed  how  to  optimize  the  size  of  a  body  with  external  forced 
convection  when  the  heat  transfer  duty  is  specified,  and  when  the  total  entropy  generation  rate  or 
the  total  cost  must  be  minimized.  We  illustrated  the  two  methods  by  modelling  the  body  as  a 
cylinder,  sphere  or  flat  plate. 

The  main  conclusion  made  visible  by  the  charts  of  Figs.  5.1  and  5.3  is  that  the  optimal 
size  increases  monotonically  with  the  heat  transfer  duty  parameter.  The  chief  objective  of  the 
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Figure  5.4  The  optimal  sizes  of  Fig.  5.3  correlated  in  terms  of  a  universal  duty  parameter  (F, 
F*)  and  the  Reynolds  number  based  on  the  transversal  length  scale  of  the  flow. 


analysis  was  to  reveal  the  most  appropriate  way  to  nondimensionalize  the  heat  transfer  duty 
parameter,  so  that  (1)  the  Prandtl  number  does  not  appear  as  an  additional  parameter  in  the  final 
charts  (Figs.  5.1  and  5.3),  and  (2)  all  the  curves  show  the  same  trend.  This  common  trend  is  an 
indication  that  a  body  with  a  more  complicated  shape  (for  example,  a  short  cylinder,  or  a 
spheroid)  can  be  optimized  numerically  based  on  a  procedure  that  will  have  the  same  steps  as 
that  presented  in  section  5.2  or  section  5.3. 

Progress  in  the  direction  of  applying  the  present  results  directly  to  the  optimal  sizing  of 
other  shapes  was  made  in  the  correlations  presented  in  Figs.  5.2  and  5.4.  A  universal  relation 
between  optimal  size  and  heat  transfer  duty  emerges  when  (1)  the  size  parameter  is 
nondimensionalized  as  a  Reynolds  number  based  on  the  transversal  length  scale  of  the  flow  (Res, 
Ren),  and  (2)  the  duty  parameter  is  nondimensionlized  as  the  total  heat  transfer  rate  of  an 
equivalent  body  in  two-dimensional  flow.  This  unique,  universal  trend  is  reported  in  Figs.  5.2 
and  5.4. 

5.5  Notation 

B 
Bs 

Ce 
Cm 
C 
Cd 
D 
F 

Fd 

Fs 

F:, 

h 
k 
L 
Nu 
Pr 

q 
q" 

Rcdjl.s 

‘-’gen 


duty  parameter,  equation  (5.3) 
sphere  duty  parameter,  equation  (5.5) 
sphere  duty  parameter,  equation  (5.8) 

unit  cost  of  electric  energy,  $/J 
unit  cost  of  material,  $/m3 
total  cost,  $ 
drag  coefficient 
diameter 

duty  parameter,  equations  (5.11)  and  (5.15) 
drag  force 

sphere  duty  parameter,  equation  (5.13) 
sphere  duty  parameter,  equation  (5.18) 

average  heat  transfer  coefficient 
fluid  thermal  conductivity 
swept  length  of  plate 
overall  Nusselt  number 
Prandtl  number 
total  heat  transfer  rate 
heat  flux 

Reynolds  number,  (D,  L,  5)  Uoo/v) 
entropy  generation  rate 
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plate  thickness 

time  of  operation 

average  surface  temperature 

free  stream  temperature 

free  stream  velocity 

plate  width,  cylinder  length 

boundary  layer  thickness,  flat  plate 

viscosity 

kinematic  viscosity 
density 

average  wall  shear  stress 
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